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Preface 


This  thesis  investigates  the  application  of  a discrete  C-Star  (C*) 
transient  response  controller  to  the  unstable  longitudinal  dynamics 
of  the  YF-16  Lightweight  Fighter  Prototype  Aircraft.  A reduced  state 
model  of  the  aircraft  is  developed  from  wind  tunnel  data  and  analyzed 
for  open  loop  stability.  This  model  is  transformed  into  a discrete 
time  domain  state  model  and  a discrete  cost  function  applied  to  de- 
velop a controller  capable  of  tracking  a commanded  response  with  zero 
steady  state  error  within  the  confines  of  a defined  C envelope  bound- 
ary. The  effects  of  both  a Zero  Order  and  First  Order  Hold  on  the 
system  C response  are  ainalyzed  using  a digital  computer  simulation. 

Topics  such  as  sample  rates,  weighting  parameters,  saturation  effects, 
and  migration  of  closed  loop  system  roots  are  also  presented. 
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Abstract 

The  YF-16  fighter  aircraft  represents  a radical  departure  from 
conventional  aircraft  design.  Reduced  longitudinal  static  stability 
resiilts  in  an  aircraft  idiich  is  unstable  in  subsonic  flight;  a 
characteristic  of  considerable  challenge  in  its  control  aspects. 

The  present  analog,  fly-by-wire  configuration  of  the  aircraft's  control 
system  makes  it  an  attractive  candidate  for  digital  control  adaptation. 

Such  a scheme,  if  successful,  could  mean  a more  compact,  lighter,  less 
fail\ire  prone,  and  more  adaptable  control  system. 

This  thesis  investigates  the  feasibility  of  a discrete  digital 
flight  controller  for  the  YF-16  through  the  design  and  analysis  of  an  I 

optimal  discrete  controller  at  M = .8  at  sea  level.  The  investigation  j 

is  limited  to  the  longitudinal  pitch  axis  only.  A reduced  state,  short 
period  approximation  mathematical  YF-16  model  is  developed  from  avail- 
able data.  The  open  loop  stability  and  response  characteristics  of  the.  1 

model  are  shown  to  be  iinacceptable , necessitating  the  use  of  closed 
loop  compensation.  The  minimization  of  a discrete  cost  function  is 
used  to  develop  a recursive  discrete  control  formula  which  uses  present 
values  of  output  and  past  error  information  to  successfully  control  i 

this  intentionally  unstable  aircraft  system.  The  thesis  discusses  ard 
incorporates  the  concept  of  a proposed  C (C-Star)  handling  qualities 
criterion  in  the  determination  of  acceptable  response.  The  concept 
of  C , which  blends  system  states,  is  incorporated  as  aui  integral  part 
of  a closed  loop,  discrete  control  model  for  the  YF-16.  Digiteil  com-  i 

I 

puter  simulation,  using  a Zero  Order  Hold  (ZOH)  or  First  Order  Hold  ! : 

(FCH)  control  scheme,  results  in  a stabilized  system  model  whose  output 


xiv 


falls  within  the  bounds  of  a defined  C*  envelope,  and  capable  of  per- 
forming the  limited  tracking  task  of  following  a 1-G  climb,  pilot, 
input  command.  T^rpical  results  in  the  form  of  plotted  time  histor7 
information  are  discussed.  Results  of  the  simulation  show  the  ZOH 
superior  to  the  FOH  control  scheme  in  reducing  the  elapsed  time  to 
reach  steady  state.  The  time  reqxiired  to  achieve  steady  state  is 
also  shown  to  be  appreciably  uneffected  at  sample  rates  greater  than 
T = 1/50  second.  More  frequent  sampling,  however,  does  result  in  the 
production  of  smaller  controls.  The  variation  in  optimal  feedforward 
and  feedback  gain  for  various  combinations  of  sample  rate  and  cost 
function  penalty  parameters  is  also  presented  in  tabular  form.  An 
increase  in  the  control  penalty  weighting  (R),  while  holding  the 
sample  rate  constant,  is  seen  as  having  the  same  effect  on  the  elapsed 
time  for  the  transient  response  to  reach  a zero  error  steady  state 
value,  as  decreasing  the  sample  rate,  while  R is  held  constant.  In 
either  case,  the  elapsed  time  increases. 

The  possibility  of  actixator  satiu^ation  due  to  excessive  control 
movement  rates  is  pointed  out  and  discussed  in  relation  to  the  simula- 
tion conducted.  Finally,  an  investigation  of  the  closed  loop  system 
complex  conj\igate  root  migrations  show  R as  being  inversely  proportional 
to  the  system  natural  frequency  ( Wn  ) , with  the  trajectory  error  penalty 
weighting  (Q)  determining  the  damping  ratio  ( { ) • 
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INVESTIGATION  OF  A DISCRETE  C-STAR 
TRANSIENT  RESPONSE  CCXJTROLLER  FOR 
THE  IF-16  AT  A SELECTED  FLIGHT  COTOITION 

I.  Introduction 

The  introduction  of  flight  control  systems  totally  based  upon  an 
electrical  primary  flight  control  system,  emphasizing  feedback,  such 
that  vehicle  motion  is  the  controlled  parameter,  is  a recent  occurrence 
in  aircraft  flight  control  system  design.  This  is  in  part  due  to  the 
infancy  of  its  technology  and  the  sense  of  security  attached  to  proven 
mechanical  flight  control  systems.  In  an  effort  to  dispel  this  resis- 
tance, and  add  to  the  literature,  but  moreover,  to  analyze  some  of  the 
problems  encountered  in  such  an  approach,  this  investigation  develops 
a linear  time-invariant  model  of  the  YF-16  Lightweight  Fighter  Proto- 
type and  combines  with  it  a digital  flight  control  system  capable  of 
tracking  step  inputs. 

Background 

The  important  credentials  of  any  would  be  air-superiority  fighter 
are  its  speed,  maneuverability  and  loiter  capability. 

Speed  is  enhanced  by  lightweight  aircraft  components  which  con- 
tribute to  favorable  thrust  to  weight  ratios;  maneuverability  is 
enhanced  by  the  aerodynamic  design  including  the  displacement  control 
surface  technique  used;  and  loiter  capability  is  enhanced  by  low  fuel 
consTimption  aind  low  failure  rates.  In  many  instances,  one  attribute 
is  enhanced  at  the  expense  of  another.  More  often  than  not,  the 

1 


r 

I 


resulting  design  reflects  a compromise  of  these  important  attributes 
and  a corresponding  compromise  in  the  resulting  aircraft  performance, 

1 VB.th  these  considerations  in  mind,  the  design  of  fighter  aircraft 

i 

becomes  a challenging  task  and  is  a subject  of  intense  interest  in  the 
Air  Force.  The  YF-16  Lightweight  Fighter  Prototype  was  designed  to 
meet  this  challenge.  It  represents  a radical  departxire  from  conven- 
tional aircraft  design  in  an  attempt  to  achieve  an  optimal  blend  of 
speed,  maneuverability,  and  loiter  capability. 

I VB.th  a length  of  just  over  forty-six  feet,  a wing-span  of  thirty- 

' one  feet,  and  a maximum  weight  of  27,000  pounds,  the  YF-16  is  cons'.der- 

I 

I ably  lighter  and  quite  a bit  smaller  than  most  present-day  fighters. 

I 

j Less  discernable  to  the  eye,  however,  the  YF-16  exhibits  a novel  center 

of  gravity  (eg),  aerodynamic  center  (ac)  relationship.  A short  dis- 
cussion of  this  important  relationship  is  in  order. 

If  the  aircraft  center  of  gravity  is  located  at  the  aerodynamic 
center,  there  is  a condition  of  zero  static  margin  or  neutral  stability 
(Fig.  l).  This  situation,  though  not  unstable,  co\ild  present  a diffi- 
CTilt  dynamic  sitviation  for  a pilot  to  control.  If  the  eg  is  moved  aft  • 
of  the  ac,  the  aircraft  becomes  unstable  (Fig.  2).  The  aircraft  will 
not  maintain  a trim  condition  auid  if  disturbed,  will  continue  to  pitch 
up  or  down  at  the  induced  rate.  The  amovint  the  eg  is  aft  or  forward 
of  the  ac  is  the  margin  of  stability.  If  the  eg  is  aft  of  the  ac,  the 
static  margin  is  termed  negative  vdiile  a eg  ahead  of  the  ac  describes 
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a.c 


Fig,  1.  Zero  Static  Stability. 


a.c. 


Fig.  2.  Negative  Static  Stability  (Unstable). 

Traditionally,  aircraft  have  been  designed  to  be  aerodynamically 
stable  (positive  static  margin)  in  the  longitudinal  mode.  In  addition 
to  allowing  for  a trim  condition,  this  traditional  design  tends  to 
bring  the  aircraft  back  to  a straight  and  level  attitude  vdien  disturbed 


from  a trimmed  flight  condition.  Were  it  not  for  this  intentionally 
designed-in  characteristic  of  static  stability,  a pilot  would  be 
continuously  adjusting  and  compensating  to  keep  the  aircraft  flying. 

In  subsonic  flight  this  type  of  stability  (positive  static 
stability)  is  desirable;  however,  as  the  aircraft  passes  through  the 
transonic  region  (Mach  .S’*"  -•  1.2”)  and  becomes  supersonic,  the  ac  moves 
aft.  Here,  its  statically  stable  nature  becomes  detrimental.  In  ef- 
fect, since  the  ac  is  now  further  aft  of  the  eg,  the  statically  stable 
aircraft  becomes  "super  stable."  This  means  that  the  aircraft  displays 
greater  resistance  to  any  disturbance  which  might  cause  it  to  move  from 
its  stable  attitude.  VJhen  such  disturbances  result  from  control  com- 
mands transmitted  by  a pilot,  the  aircraft,  not  being  able  to  discrim- 
inate between  a gust  or  legitimate  command,  resists.  This  is  an 
obvious  handicap  for  a modem  high-performance  fighter  which  is  expected 
to  display  exceptional  agility  in  maneuvers  at  supersonic  speeds. 

The  obvious  solution  to  this  problem  is  to  reduce  the  longitudinal 
stability  margin  in  order  to  make  the  aircraft  more  maneuverable  at 
supersonic  speeds.  Such  a solution,  however,  would  cause  an  aircraft 
to  be  tanstable  in  subsonic  flight.  Such  a situation  could  not  be 
tolerated.  This  obvious  dilemma  existed  for  sometime,  until  recent 
developtoents  in  automatic  control  system  technology  offered  a solution. 

The  solution  to  this  dilemma  takes  the  fom  of  what  is  called  a 
fly-by-wire  (FBW)  control  system.  Instead  of  a complex  network  of 
mechanical  linkages,  this  system  uses  electronic  signals  to  relay 
requested  response  requirements  from  the  pilot  to  electro-hydraiolic 
servos  which  move  the  control  surfaces.  In  a constant-G  climb,  for 
instance,  as  the  aircraft  response  to  a command  begins,  the  response 
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is  fed  back  to  the  flight  control  computer  where  it  is  compared  to  the 
pilot's  G requirements.  TOien  the  two  match,  the  signal  is  milled;  no 
further  control  is  applied.  The  aircraft  maintains  the  constant  G 
climb  until  the  command  is  changed. 

The  IF-16  uses  such  a system.  However,  the  relaxation  of  static 
stability  and  all  its  performance  benefits  have  not  been  without  cost. 

Since  the  IF-lb  is  intentionally  designed  to  maintain  from  seven  to  ten 
percent  negative  static  margin,  the  stability  of  the  aircraft  is  reduced 
to  the  extent  that  the  YF-16  is  unstable  in  pitching  motion  at  all  sub- 
sonic airspeeds.  Therefore,  with  an  airframe  intentionally  designed  to 
be  \instable  in  subsonic  flight,  the  IF-16  relies  upon  the  flight  con- 
trol system  to  maintain  tight  closed-loop  control  over  its  unstable 
dynamics. 

The  fly-by-wire  control  system  is  used  to  keep  a tight  rein  on  the 
aircraft,  especially  during  the  subsonic  portion  of  the  flight  envelope 
>diere  it  is  inherently  unstable.  Because  of  the  implementation  of  re- 
laxed static  stability  technology,  the  flight  control  systems  acts  as 
a compensator  to  fill  in  for  the  inadequacy  of  the  aircrad!t  design  in 
terms  of  stability.  The  end  result  is  that  the  IF-16  is  not  only 
stabilized  but  has  the  advantages  of  reduced  drag  as  a result  of  the 
unstable  airframe  design  and  greater  ease  of  maneuverability  especially 
at  supersonic  speeds.  As  the  gains  and  compensation  networks  in  the 
control  system  are  changed,  the  stability  of  the  aircraft  can  be  varied. 

The  advantage,  then,  is  that  in  subsonic  flight,  the  adrcraft  can  be  | 

made  stable,  to  fly  in  the  conventional  manner,  and  in  supersonic 
flight,  the  stability  can  be  relaxed  to  allow  greater  maneuverability. 
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At  present,  the  controller  for  the  unstable  pitch  axis  is  a complex 
analog  computer.  Gain  adjustments  are  controlled  by  this  on-board  de- 
vice, Gains  are  "scheduled"  as  a function  of  flight  condition  in  an 
adaptive  control  scheme.  For  exan^Jle,  one  particular  airspeed, 
altitude,  angle  of  attack  and  pitch  rate  might  identify  a particular 
gain  from  the  schedule  of  possible  gains,  while  a slightly  different 
airspeed,  angle  of  attack,  and  pitch  rate  might  identify  another.  In 
its  eUialog  configiiration,  the  present  controller  is  not  only  bulky  and 
difficult  to  mechanize,  but  also  presents  problems  in  terms  of  future 
modification. 

With  some  appreciation  for  the  present  control  configuration  of 
the  XF-16  just  described,  it  has  been  suggested  that  digital  computers 
be  substituted  for  present  analog  computer  controllers  (Ref,  5).  This 
3\;ggestion  is  reasonable  in  light  of  the  advent  of  small  digital  com- 
puters and  recent  advances  in  digital  integrated  circuits  which  make 
the  digital  computer  an  attractive  candidate  as  a controller.  Much 
effort  has  been  expended  in  researching  this  proposed  substitution 
approach  (Ref,  18), 

A digital  controller  would  have  the  advantage  of  being  more  com- 
pact, of  less  mechanical  complexity,  lighter,  less  failure  prone,  and 
more  adaptable  to  changes  in  the  form  of  software  adjustments  to  modify 
control  laws.  However,  with  these  advantages,  are  associated  new  prob- 
lems not  encoiaitered  with  the  analog  controller.  Two  such  problems, 
among  a list  of  others  including  the  areas  of  sufficient  word  length 
and  adequate  memory  size,  are  the  identification  of  an  adequate  recvir- 
sive  digital  control  algorithm  and  the  determination  of  the  sampling 
rate  with  which  to  implement  this  algorithm. 
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Problem  Statement 

Along  this  vein,  the  problem  addressed  in  this  study  is  the  design 
and  analysis  of  an  optimal  discrete  controller  at  a selected  flight  con-  i 

dition  for  the  longitudinal  pitch  axis  of  the  YF-16  using  a specific  C* 

(pronounced  C-Star)  performance  criteria.  The  scope  of  the  analysis  is 
limited  to  the  consideration  of  the  unstable  longitudinal  mode  of  the 
TF-16  only.  It  is  this  mode  which  is  applicable  to  the  C performance 
criteria.  Additionally,  it  is  pitch  response  which  is  the  most 
divergent  dynamic  mode. 

The  cases  of  both  Mach  .8  and  Mach  1.2  at  sea  level  and  30,000 
feet  for  trim  angles  of  attack  are  considered.  These  particular 
flight  conditions  were  selected  based  on  the  availability  of  data 
(Ref.  13  and  14),  and  in  consideration  of  the  fact  that  high  subsonic 
flight  at  sea  level  is  considered  the  most  critical  area  for  pitch 
response  for  this  aircraft.  Here,  control  surface  deflections  produce 
the  greatest  dynamic  effect  on  the  statically  unstable  airframe. 

Order  of  Presentation 

In  an  effort  to  develop,  simulate,  and  discuss  a suitable  discrete 
control  concept  for  a pitching  model  of  the  YF-16,  the  remaining  chap- 
ters of  this  investigation  are  organized  as  follows. 

In  Chapter  II,  non-dimensional  stability  derivatives  for  Mach  ,8 
at  sea  level  for  a three  degree  of  freedom  aircraft  model  are  developed. 

The  corresponding  derivatives  for  the  three  remainiTig  flight  conditions 
are  also  listed.  The  model  is  derived  from  available  wind  tunnel  data. 

In  Chapter  III,  the  equations  of  motion  developed  in  the  preceeding 
chapter  are  investigated.  The  nature  of  the  open  loop  system  stability 
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and  the  character  of  its  response  performance  are  discussed. 


i 

! In  Chapter  IV,  a discrete  model  of  the  proposed  control  system  is 

1 

presented.  Additionally,  the  concept  of  C as  a plausible  stability 
I ciateria  is  presented.  This  approach,  which  has  gained  some  popularity, 

uses  a linear  blend  of  normal  acceleration,  pitch  rate,  and  pitch  accel- 
eration to  define  proper  performance.  The  C*  concept  is  included  as  an 
integral  part  of  the  proposed  control  system. 

Chapter  V details  the  development  of  a discrete  system  model  from 
its  continuous  representation.  Also  included  is  the  development  of  the 
observer  matrix  needed  to  transform  the  system  states  into  usable  C 
parameters. 

Chapter  VI  presents  the  mathematical  development  of  the  optimal 
discrete  controller  using  this  C performance  criterion  and  the  rec\ir- 
sive  digital  control  algorithm  which  results.  Also  discussed  is  the 
concept  of  a hold  device.  The  chapter  concludes  with  a discussion  of 
the  structvire  of  the  software  developed  for  this  investigation  and  i 

used  to  implement  the  controller  concept.  This  program,  as  a function 
of  sample  rate,  determines  the  optimal  gains  to  be  applied  in  the  con- 
trol algorithm  to  achieve  satisfactory  transient  C*  response  of  the 
system.  Various  assumptions  and  limitations  are  also  discussed  along 
with  the  simulation  technique  conducted.  The  controller  is  simulated 
at  various  sample  rates  using  a CDC  6600  Digital  Computer.  The  migra- 
tion of  the  closed  loop  system  roots  as  a function  of  sample  rate  is 
investigated  along  with  saturation  effects  on  the  horizontal  stabilizer 
servo.  Using  the  simulation,  an  attempt  is  made  to  identify  a minimum 
sample  rate  which,  based  upon  the  satiiration  limit  of  the  pitch  control 
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servo  and  characteristic  root  locations,  still  results  in  a response 
within  defined  C envelope  bounds. 

Chapter  VII  discusses  observations  on  the  controllability  of  the 
system  resulting  from  the  simulation  effort,  lypical  results  are 
presented  and  explained.  The  effects  of  sample  rates  and  cost  function 
penalty  weighting  variations  are  also  included. 

Finally,  in  Chapter  VIII,  the  conclusions  drawn  from  this  inves- 
tigative effort  are  presented  and  recommendations  are  discussed 
pertaining  to  fiirther  research. 
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H.  Mathematical  Model  and  the  Reducation  of  Data 

« 

Before  the  implementation  of  any  controller  and  an  analysis  of 
its  effectiveness  can  be  undertaken,  an  adequate  mathematical  model 
of  the  aircraft  must  first  be  developed.  The  purpose  of  this  chapter 
is,  therefore,  to  develop  a mathematical  model  of  the  YF-lS  for  longi- 
tudinal pitching  motion.  Generalized  equations  of  motion  are  first 
presented  and  then  reduced  in  complexity  through  the  use  of  simplifying 
assumptions.  The  remainder  of  the  chapter  presents  the  development  of 
rigid  stability  derivatives  which  constitute  the  individual  terms  of 
the  equations  of  motion.  These  derivatives  are  summarized  in  table 
format  at  the  conclusion  of  the  chapter. 

Eqmtions  of  Motion 

As  explained  in  the  preceding  chapter,  this  investigation  is 
concerned  only  with  the  longitudinal  pitching  mode  of  the  IF-16. 
Therefore,  disregarding  the  negligible  cross-coupling  effects  of 
lateral-directional  motion,  the  small  perturbation  equations  of  motion 
which  describe  the  dynamics  of  the  aircrsift  follow: 
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= c. 


(3) 


These  equations  are  Blakelock's  non-dimensional  equations  of  motion 
for  the  longitudinal  axis  of  an  aircraft  and  vriLU  be  used  to  model 
the  YF-16.  Their  development  vdll  not  be  presented  here  but  can  be 
found  in  Reference  1.  Definitions  of  individual  equation  elements 
are  included  in  the  previous  list  of  Abbreviations  and  Symbols,  p.viii. 

The  various  non-dimensional  coefficients  in  these  equations  are 
referred  to  as  stability  derivatives.  The  equations  can  be  simplified 
some\diat  by  the  judicious  elimination  of  three  of  these  derivatives. 

It  is  possible  to  eliminate  > and  terms  from  the 

equations  thus  considerably  simplifying  the  modeling  equations. 

As  Blakelock  points  out,  and  , the  effect  of  downwash 
from  the  vring  on  the  horizontal  tail  and  the  effect  of  pitch  rate  on 
drag,  respectively,  have  negligible  contributions  to  the  equations 
and  are  usually  neglected. 

, a term  resulting  from  slipstream,  thrust,  and  flexibility 
effects,  is  the  change  in  pitching  moment  due  to  a change  in  forward 
velocity.  As  pointed  out  in  Reference  1,  this  term  can  be  safely 
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neglected  for  jet  aircraft.  Using  these  simplifications  for  all  the 
flight  conditions  in  this  sttidy  results  in  equations  (l),  (2),  and 
(3)  being  expressed  as; 


- 

• 

• 

ksi-c  c. 

h - 

«c 

+ 

<S>J  e 

J 

- 

(4) 


O 

A‘k 


(6) 


As  Reference  1 points  out,  the  forcing  function  of  Cp  can  be 

*a 

approximated  as  equaling  zero,  while  Cp  and  CV  equal 

ie  , respectively.  Now,  substituting  for  the  forcing  ibnctions 
on  the  right-hand  side  of  equations  (4),  (5),  and  (6),  with  the  nota- 
tional  difference  that  , the  deflection  of  the  horizontal  stabilizer, 
replaces  , the  deflection  of  the  elevator,  since  the  IF-16  .has  a 
movable  horizontal  stabilizer  and  no  elevator,  the  equations  become: 
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In  order  to  evalviate  and  apply  these  equations  to  the  YF-16,  the 
various  terms  must  be  replaced  by  actual  values.  The  remainder  of  this 
chapter  presents  the  detailed  development  of  each  of  these  stability- 
derivatives  for  the  flight  conditions  at  Mach  .8  at  sea  level.  Flexible 
aircraft  effects  will  not  be  considered. 

As  a prelude  to  this  development,  the  following  directional 
convention  is  presented. 

Directional  Convention 

A positive  control  force  produces  a negative  siirface  deflection 
and  causes  a positive  moment  about  the  pitch  axis  (Fig.  3). 


Note  that  for  the  horizontal  stabilizer  control  s-urface,  trailing 
edge  up  is  defined  as  negative,  i.ei,  ^ pitch  up. 

The  rigid  stability  derivative  calculations  are  based  on  stability 
and  flight  control  data  furnished  by  General  Dynamics  Convair  Aero- 
space Di-vision  (Ref.  13  and  14). 
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The  configoiration  chosen  was  that  of  a clean  aircraft  (i.e., 
gear  up,  flaps  retracted,  etc.)  with  external  fuel  tanks  removed. 
The  physical  specifications  and  atmospheric  conditions  listed  in 


Table  I 

Physical  Specifications 

Nominal  Weight  (YTT) 

16,519  lbs 

Wing  Area  (S) 

280  ft^ 

Mean  Aerodynamics  chord  (c) 

10.937  ft 

Mass  Moment  of  Inertia  (lyy) 
Distance  from  c/k  of  wing  MAC 

V 

39,199  ft-lbs-sec^ 

to  c/4  of  Horizontal  Tail  MAC 

15.66  ft 
(See  Appendix  A) 

Table  II 

Atmospheric  Conditions 

Altitude 

Sea  Level 

30,000'  (FL  300) 

Mach  No. 

.8 

1.2 

.8  1.2 

1 

Density  ( p ) o 
slu^ff^ 

.002378 

.002378 

.00890  .000890 

Dynamic  Pressiire  (q) 
lb/ft2 

949.44 

2136.24 

282.02  634.51 

Velocity  ( ^ ) 
ft/sec 

893.6 

1340.4 

796.08  U94.1 

As  a preliminary  step,  it  was  necessary  to  determine  the  trim 
angle  of  attack  for  the  selected  flight  condition  (M  = ,8,  sea  level). 

Trim  Angle  of  Attack 

At  the  trim  condition  #g's  = 1.0, 


C^9o) 


I.O 


/.o 


» .73* 


(13) 


(14) 


Here,  it  has  been  assumed  that  the  pitching  moment  and  drag  equations 
are  satisfied  at  this  trim  angle  of  attack.  This  angle  of  attack  was 
used  to  extract  data  from  the  various  graphs  and  tables  found  in 
References  13  and  14. 

Rigid  Stability  Derivatives 

The  purpose  of  this  section  is  the  determination  of  exact  values 
of  the  stability  derivatives  at  M = .8  at  sea  level  needed  to  mathe- 
matically model  the  aircraft.  These  derivatives  are  based  on  a 
stability  aixis  system  with  its  centroid  at  the  wing  quarter  chord 
and  waterline  of  91.0  inches.  Detailed  definitions  of  each  derivative 
beyond  the  short  definitions  found  in  the  List  of  Symbols  on  page  viii, 
are  found  in  Reference  1 and  will  not  be  included  here. 
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M = 533.45  slugs 


Forward  Vel  (^): 

Speed  of  sound  at  sea  level 

c 

= 1117  ft/sec 

(18) 

Mach  .8 

= .8  (1117) 

(19) 

'U 

= 893.6  ft/sec 

17 


e.  Dynamic  Pressure  (q); 


q - — £ Jt*- 

A 

= I.002378  (893.6)^ 


q = 949. A4  Ib/ft^ 

f iJL  . 

hU  ^ (513. 45)  (893. 6) 

Sj,  (280)(949.44) 


= 1.726 


39,199 

280  (849.44) (10.937) 


h 


_i . 

av 

s. 

Qli 


e 


i. 


Cw 


.0135  sec 


2 


10.937 
2 (893.6) 


.0061  sec 


“Mg  = -513.45  (32.1725) 

Sq  (280)  (949.44) 


= -.0621 


(20) 


(21) 


(22) 


(23) 


(24) 
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was  determined  for  pitch  up  only.  Data  was  available 
for  horizontal  tail  deflections  of  0®  -10“  which,  based 

upon  the  previously  set  convention,  causes  the  aircraft's 
nose  to  pitch  up. 


(25: 


A linear  interpolation  of  the  data  between  od  = 0“  aind 
= 2.5®  was  performed  to  determine  Cj^  for  = .73 


06 

0 

0 

II 

.1120 

II 

0 

• 

.1126 

2.50° 

.IIW) 

1126 


(57.3  deg) 


This  viG.ue  was  extracted  directly  from  the  tables. 


C = ^.3900 


1. 


-i— 

it  i 

-3.0647 


m. 


.6981  • (-4.3900) 


- C 


(28) 


(29) 


(30) 


for  = .73  C.  = .OOn/deg  = .0630/Rad 

T < 

= .0625  (interpolation 
resiilt) 


= .0625  - .0630 

= -.0005 


(31) 

(32) 

(30) 


2r  - ic. . a 


(33) 
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The  first  term  on  the  right-hand  side  of  equation  (33) 

9t 

reduces  to  zero,  that  is  — = 0,  since  thrust  is 
essentially  constant  for  jet  aircraft.  Equation  (33) 
can  now  be  expressed  as; 


c,  = - iC^  - U 

or  equivalently  as: 


c,  = - a Ck 


.02Z^8  (interpolation  result) 


= .0496 


= f^\  .0248  - .0268 

* L -8  - -9 


= .0160 


= -.0656 


= -.0496  - .0160 


o.  : 


C - 


c,  . 
« 


- " ^1.. 


- = -0.0248 


<^1.  = (57.3)(.078) 


= 4.4694/Rad 


C - - c - 0 


= -.0248  - 4.4694 


= -4.4942 


P-  : 


= J_  C 
# ^**»c 

% e 


= (.6981) (-.6452) 


-.4504 
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If  positive,  the  aircraft  is  statically  mstable. 

Recalling  the  discussion  of  stability  in  the  previous  chapter, 
we  woiHd  expect  the  YF-16  to  have  a positive  for  this 

subsonic  flight  condition.  As  shown  earlier; 


= 4.4694/Rad 

(42) 

^ " ^.c.^ 

(47) 

= ^ .35  - .3^ 

S.M.  = .04 

(45) 

= (.04)(4.4694) 


= .1788 

d 
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As  anticipated,  = .8,  S.L.  is  positive.  A check  of 

n 

M ~ 1.2,  S.L.  shows  this  coefficient  as  negative. 


= (S.M.) 

M*<-»  Sa. 


(48) 


= (.35  - .56)(5.2143) 

= -1.0950 

•t 

t.L. 

For  this  stability  derivative  then,  the  model  truly  reflects  the 
variable  stability  aspect  of  the  YF-16,  i.e.,  a shift  from  in- 
stability to  a stable  nature  as  the  aircraft  accelerates  past 
M = 1.0. 


c. 


Cy  = - Ol  C 


(49) 


or  equivalently  as: 


aC,, 

Cj  = - - NacI'  A^eit 

^ L • 


(50) 


As  earlier. 


Cl  = .0625 


and 


-AC^=  -.1250 


h 


aC,. 


Ah 


^ M..g  * ^4.  ^/VC.9  j 


(32) 

(51) 

(52) 
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41 


I 


i 

i 


t 

f 


= -,1250  - .0632 

Cr  = -.1882 

A tabialation  of  the  values  just  calculated  to  be  used  in  equations 
(10),  (11),  nd  (12)  follows  in  Table  III.  Also  involved  are  the  sta- 
bility derivatives  for  the  three  remaining  flight  conditions.  These 
values  were  calculated  in  the  same  manner  as  the  case  of  M = .8  at 
sea  level,  just  presented. 

Using  the  appropriate  data  from  Table  III,  equations  (lO),  (ll), 
and  (12)  can  be  expressed  numerically. 


YF-16  Longitudinal  Equations  of  Motion 


M = .8,  Sea  Level 


|1.726s  + .065^  U(s)  + .0005^  ^ (s)  + 

'!o62i]  ©(s)  = 0 

Q .1882  ^i(s)  +Q.7325S  + 4.494^  *<(3)  - 

[1.70733  ®(s)  = .4504  j^(S) 


(53) 


(54) 
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|To093s  - .1788j<i<(s)  + |^.0135s^  + .0268sj  0(s)  = 


M = 1,2,  Sea  Level 

[l.l506s  + .103^  ^i(s)  + .0813  J'<(s)  + [^.0276^  0(s)  = 0 


(56) 


1 ^ r 

.1200 

(s)  + ] 

J ^ 

^.8132  JoC 


•^s)  - 


[l.ll862s]  ©(s)  = -.36iS'^(s) 


(57) 


|^-.00l64s  + .7392jii(s)  + .006s2  + ,oi855s  ©(s)  = 


Analysis  of  these  equation  sets  is  the  subject  of  Chapter  III. 

Investigation  of  the  corresponding  sets  of  equations  for  the 
remaining  two  flight  conditions  at  30,000  feet  shows  a close  corre- 
spondence in  root  locations  to  the  roots  of  the  sets  of  equations 
above. , Due  to  this  situation,  the  remainder  of  the  investigation 
vtn  concentrate  on  the  above  sets  of  equations  leaving  the  higher 
altitude  flight  cases  for  reference  and  futtore  investigation. 


III.  Analysis  of  the  Eqiiations  of  Motion  and  System  Stability 


Having  derived  two  sets  of  equations  of  motion  in  the  preceding 
chapter,  the  IF-16  model  can  now  be  analyzed  regarding  its  stability 
and  response  performance.  To  accomplish  this,  this  chapter  begins 
with  the  determination  of  the  characteristic  equations  for  flight  at 
M = .8,  and  1.2  at  sea  level.  For  a given  system,  there  is  only  one 
polynomial,  termed  the  system  characteristic  equation,  which  determines 
the  form  of  the  system  transient  response  regardless  of  the  type  of 
signal  chosen  as  the  input.  This  polynomial  is  determined  and  its 
roots  investigated.  Short  period  approximations  of  both  flight  con- 
ditions, vdiich  assume  zero  pertxu’bation  in  forward  velocity  as  the 
aircraft  maneuvers,  are  then  developed  and  the  roots  of  the  associated 
characteristic  polynomials  investigated.  Open  loop  system  responses 
to  various  inputs  are  discussed.  After  a brief  discussion  on  the  high 
frequency  nature  of  the  aircraft  actuator  servo,  the  chapter  concludes' 
with  a variable  gain  root  locus  analysis  of  system  stability  for  both 
flight  conditions. 

For  reference,  the  modeling  eqioations  are  repeated  here. 


M = .8.  Sea  Level  (Set  A) 

Jl.726s  + .0656]  <l(s)  +[  .0005]«<(s)  + [.062lJ  © (s)  = 0 (53) 

J.I882]  U(s)  + [ 1.7325s  + 4.49A2]«<(s)  - [l.7073s]  ©(s)  = 

-.4504  (54) 


J.OO933  - .1788]i(s)  + ^ .01358^  + .0268aJ  ©(s)  = 


-.6452  S^(s) 


(55) 
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M = Sea  Level  (Set  B) 


^1.1506s  + .103^U(3)  .081^ti(s)  + 

.0276]  ©(s)  = 

0 

(56) 

|^.12ju(s)  + ^1.14653  + 4.8I32]  «ic(s)  - 

1.11862s]  ©(s)  = 

-.361 

(57) 

[*-.001643  + .7392lot(s)  + [ .006s^  + .01855s  0(s) 

W — < — 

-.5157 

(58) 

The  first  set  of  three  equations  will  be  referred  to  as  Set  A while 
the  latter  as  Set  B. 

AnaTysls  of  Equations  of  Motion 

The  characteristic  equations  for  both  sets  of  equations  were 
determined  using  a digital  computer  subroutine  specifically  designed 
for  control  systems  analysis  (Ref.  l6).  Details  of  the  formatted  in- 
structions and  computer  results  are  presented  in  Appendix  B.  The 
results  of  this  analysis,  the  model  characteristic  equations  auid 
respective  roots,  are  as  follows: 

Characteristic  Equations 
Set  At 

.0403693^  + .2137993^  - .310934s^  - .012018s  - .002090  = 0 (59) 

Set  Bt 

.0079153^  + .0562993^  - 1.0590763^  + .094454s  + .002448  = 0 (60) 
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The  negative  terms  in  equation  (59)  are  indicative  of  roots  in 
the  right  half  of  the  S-plane  corresponding  to  an  unstable  situation. 
This  is  not  the  case  in  equation  (60)  where  all  roots  appear  positive. 
This  is  confirmed,  as  the  results  of  Appendix  B show,  by  factoring  the 
characteristic  equations  into  their  corresponding  roots. 


Characteristic  Roots  of  Set  A 

Real 


Tmag 


01 

0. 

01 

.7804015E 

- 01 

01 

-.7804015E 

- 01 

01 

.6967083E 

- 30 

Characteristic  Roots  of  Set  B 


Real 


..4474255E  - 01 
..4474255E  - 01 
■.3511721E  + 01 
-.351172IE  + 01 


.1790868E  - 01 
-.1790868E  - 01 
.1099289E  + 02 
-.1099289E  + 02 


These  system  root  locations  can  be  shown  pictorially  in  the  S-plane  of 


The  root  in  the  right  half  plane  confirms  the  unstable  nature  of 
the  system  described  by  Set  A.  The  conjugate  roots,  close  to  the 
imaginary  axis,  are  conventional  phugoid  roots  due  to  the  low  fre- 
quency and  correspondingly  long  period  associated  with  their  oscil- 
latory motion,  vdiile  the  two  remaining  aperiodic  roots  (Ref.  4)  are 
considered  as  short  period  roots. 


For  the  Set  B flight  condition,  unlike  Set  A,  the  roots  indicate 
that  the  system  is  stable  (Fig.  5).  Again,  the  phugoid  roots  lie 
close  to  the  imaginary  aods  while  the  remaining  two  roots  appear  as 
normal  short  period  roots.  Such  a shift  in  system  characteristics, 
from  the  unstable  case  of  Set  A to  the  stable  situation  of  Set  B,  is 
due  to  the  variable  location  relationship  of  the  center  of  gravity 
and  aerodynamic  center,  as  a function  of  Mach  number,  discussed  in 
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Chapter  I.  This  relationship  is  expressed  in  the  stability  derivative 
C|„^ discussed  earlier  and  is  the  key  term  in  determining  the  stability 
of  this  system  of  equations. 


Since  in  both  cases,  the  phugoid  roots  are  well  behaved  and  stable, 
a short  period  approxiniation  disgarding  the  phugoid  oscillatory  mode 
is  svifficient  to  adeq\iately  model  the  system.  Additionally,  it  is  the 
rapid,  transient  response  of  the  aircraft,  represented  by  the  short 
period  mode,  which  is  of  interest.  Again,  based  upon  Reference  1, 
the  short  period  equations  can  be  expressed  as  follows: 


(61) 


(62) 


The  Set  A eqiiations  reduce  to  the  following: 


(1.726s  + 4.4942)  u (s)  + (-1.726s)  0(s)  = -.4504^^(s)  (63) 

(.009272s  - .1788) «<(s)  + (.01353^  + .026779s)  0(s)  = 

-.6452jJ^(»)  (64) 

likewise,  the  Set  B equations  reduce  to: 


(1.1506s  + 4.8132)  «i(s)  + (-1.1506s)  ©(s)  = -.36l^,<i)  (65) 

(-.00164s  + .7392)  «i(s)  + (.006s^  + .018553s)  ©(s)  = 

-.51575|(«>  (66) 


The  TRANFUN  analysis  program  (Ref.  16)  is  used  in  a similar 
manner  as  presented  in  Appendix  B.  Results  show  the  characteristic 
equation  of  the  short  period  system  described  by  equations  (63)  and 
(6U)  to  be: 

.023301s^  + .1228953^  - .1882536s  =0  (( 


which  produces  the  following  roots: 


0. 

.1240223E  + 01 
-.6514491E  + 01 


In  the  S-plane,  these  appear  as  shown  in  Figure  6. 


Fig.  6.  Short  Period  Approximation  Roots  (Set  A) 


Again,  the  instability  is  present  while  the  root  locations,  due 
to  the  simplification,  have  migrated  less  than  one  percent.  The  pole 
at  S = 0 results  from  the  fact  that  0)  was  taken  as  zero  in  equations 
(61)  and  (62)  which  effectively  eliminates  the  effects  of  gravity. 

The  system  transfer  functions  of  the  reduced  Set  A equations  (63) 
and  (6U)  become: 
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-.26095  (s  + 185.132) 

(s  - 1.240223)(s  + 6.51A491)  (68) 

-47.61337  (s  + 2.686213) 

s (s  - 1.240223)(s  + 6.514491)  (69) 

The  same  process  can  be  performed  on  Set  B equations  (65)  and  (66). 

The  resulting  short  period  characteristic  equation  is: 

.0069s^  + .048343^  + .93982s  = 0 (70) 

The  roots  of  this  eqiiation  are: 

Real  Imag 

0.  0. 

-3.501  -11.130 

-3.501  +11.130 

In  the  S-plane,  these  appear  as  sho\fli  in  Figiire  7. 


M = 1.2,  S.L. 

r 

/< 

X-- 

An 

1 

-M 

1 

1 

X- 

Fig.  7.  Short  Period  Approximation  Roots  (Set  B) 
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As  expected,  the  system  remaijis  stable  with  only  a minor  shift 
in  the  root  locations.  Also,  an  additional  root  at  the  origin  is 
introduced  for  the  reason  given  earlier. 

ii 

The  two  short  period  approximation  transfer  functions  of  Set  B 
equations  (65)  and  (66): 

(3  + 277.0^7)  , , ^ 

(s  + 3.501021  + 11.13j)  (71) 


(72) 


(73) 


Response  to  Inputs 

A further  appreciation  for  the  dynamic  characteristics  of  the 
system  can  be  gained  by  looking  at  the  open-loop  system  response  to 
various  inputs. 

The  transfer  function  (Fig.  8)  is  used  to  demonstrate  this 
response.  This  particxilar  transfer  function  is  chosen  since  pitch 
(e),  and  ultimately  pitch  rate  (0),  are  of  interest  in  investigating 
the  longitudinal  dynamics  of  an  aircraft.  The  digital  computer  pro- 
gram PARTL  is  used  to  calo\ilate  the  various  responses  (Ref.  10). 
Figure  9 shows  the  pole-zero  locations  of  the  basic  airframe  for  the 
M = ,8  condition  at  sea  level  reflected  in  equation  (69). 


©(1^  _ -86.03576  (s  + 3.729762) 

S^ts)  s (s  + 3.501021  + 11.13j) 


where 


= 11.667  Rad/sec 

= 11.130  Rad/sec 

Sff.  = .3 

"^1-  = .539  sec 


become : 
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The  instability  of  the  model  in  this  flight  condition  is  evident  from 
the  pole  which  lies  in  the  right  half  plane.  Figure  10  shows  the 
pole-zero  configuration  based  on  equation  (72)  for  the  M = 1.2  case 
at  sea  level.  In  this  flight  condition,  the  stability  of  the  model 
is  reflected  by  the  fact  that  all  roots  lie  in  the  left  half  plane. 


Fig.  10.  Open  Loop  Pole  Zero  Locations 


Noti:ig  that  F(t)  represents  0 in  the  plots  which  follow.  Figure  Ua 
shows  the  response  to  a step  horizontal  stabilizer  deflection  for 
the  M = .8  at  sea  level  case;  while.  Figure  Ub  shows  the  corres- 
ponding response  for  the  high  Mach  case. 

The  need  for  some  type  of  control  is  clearly  apparent  in  both 
cases.  Figure  11a  shows  a very  rapid  parabolic  divergence  in  the 
system  for  the  M = .8  case.  Figure  11b  shows  a lightly  damped 
oscillation  superimposed  on  a ramp  type  response.  Here,  although 
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the  response  to  the  constant  step  input,  as  expected,  is  ramp- 
like  and  stable,  the  light  damping  indicated  by  eqiiation  (73) 
to  be  ^ = .3,  is  unsatisfactory.  These  response  characteris- 
tics are  more  obvious  in  Figure  12.  In  this  set  of  figures, 
the  input  is  similar  to  a unit  impulse.  This  input  is  given 
physical  significance  by  imagining  a pilot  pulling  back  on  the 
control  stick  slightly  and  then  immediately  returning  the  stick 
to  neutral  to  mil  1 the  command.  Such  a rapid  series  of  commands 
should  cause  the  nose  to  pitch-up  as  the  aircraft  attempts  to 
begin  a climb  and  then  quickly  lowering  as  the  command  is 
removed. 

Figtire  12a,  reflecting  the  low  Mach  condition,  is  very 
similar  to  Figure  Ua.  Again,  a rapid  parabolic  divergence 
results  as  the  aircraft  pitches  up  and  continues  to  do  so  even 
when  the  command  is  removed.  This  phenomenon  definitely  re- 
flects a dynamically  unstable  condition.  Figure  12a,  however, 
reflects  the  stable  nature  of  the  high  Mach  condition  since  the 
aircraft  returns  to  a stable  steady  state  after  the  perturbing 
command  is  removed.  The  extremely  light  damping  ( C = ,3)> 
more  evident  here  than  in  Figure  Ub,  remains  unsatisfactory. 

For  both  flight  conditions,  these  plots  show  the  response  of 
the  open-loop  is  imacceptable.  Some  type  of  feedback  control 
scheme  is  definitely  necessary. 
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Fig.  12a.  Imptilse  Response  at  M = .8. 

Impulse  RESPor^  rt  n-i.aTTEfl'LEvfL 




Fig.  12b.  Lnptilse  Response  at  M ~ 1.2. 
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Aircraft  Actuator  Servo 

Before  proceeding  any  further,  the  dynamics  of  the  aircraft 
servo,  which  positions  the  horizontal  stabilizer,  must  be  considered. 
Reference  13  indicates  that  the  command  servo  transfer  function  for 
this  aircraft  is: 


(52)^ 

+ 2(.7)52s  + 52^ 


while  that  of  the  power  actuator  is: 


^i..w 


s + 20 


These  are  integrated  with  the  aircraft  dynamics  as  shown  in  Figure  13. 


Command  Servo 


Power  Actuator 
Servo 


Fig.  13.  Block  Diagram  of  Servo  and  Aircraft. 


The  open  loop  transfer  fijnction  of  the  cascaded  networic  becomes: 


_ 

^ ' 


K gain  (Zero  of  A/C  Dynamics) 

[s  + 36.4  + 37.135j)(s  + 20)(Poles  of  A/C  Dynamics) 


The  conjiigate  roots  produced  by  the  command  servo  (-36.4  + 37.135j) 
are  to  the  extreme  left  in  the  S-plane.  Their  effect  on  the  system  cam 
be  safely  disregarded  because  of  their  high  frequency.  Their  only  ap- 
preciable effect  woxiLd  be  the  introduction  of  a small  phase  lag  in  the 
system.  Therefore,  the  command  servo  will  not  be  considered  fvirther 
in  the  remainder  of  this  investigation. 

Having  accomplished  this  simplification,  the  system  reduces  to 
the  closed  loop  representation  shown  in  Figure  14. 


The  effect  of  the  addition  of  the  actuator  servo  to  the  system 
can  now  be  ainalyzed.  Use  of  the  ROOTL  computer  program  (Ref.  11) 
produces  a variable  gain  root  locus  for  each  flight  condition.  The 
Mach  .8  resiilts  appear  in  Figure  15  while  Mach  1.2  results  appear 
in  Figure  l6. 
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The  results  of  this  root  locus  analysis,  in  which  a rudimentary, 
unity  feedback  system  is  implemented,  indicate  that  the  high  Mach 
case  can  be  stabilized  for  a limited  range  of  gain.  Figure  15  points 
up  a different  situation.  Here,  no  value  of  gain  will  stabilize  the 
system.  Compensation*  , beyond  more  gain  adjustment,  is  required  for 
this  flight  condition. 

Sirngnary 

This  chapter,  through  the  development  and  investigation  of 
system  characteristic  equations  and  system  input  responses  has  shown 
the  open  loop  dynamic  character  of  the  'YF-16  modeling  equations  to 
be  unsatisfactory.  A short  period  approximation  of  the  system, 
varying  only  slightly  from  a fully  developed  model  considering  both 
phugoid  and  short  period  oscillating  modes,  was  shown  to  adequately 
represent  the  system.  Additionally,  a root  locus  analysis  pointed 
out  the  need  to  compensate  the  system  especially  for  M = .8  flight 
at  sea  level. 

It  is  the  compensation  of  this  critical  flight  condition  re- 
flected in  the  Set  A equations,  vdiich  will  be  the  focus  of  attention 
for  the  remainder  of  this  investigation.  The  application  of  a classi- 
cal lag  compensator  for  analog  implementation  might  s\ifficej  however, 
as  intimated  earlier,  investigation  will  focus  on  the  implementation 
of  a discrete  optimal  controller  scheme  specifically  designed  for 
digital  computer  implementation. 

Compensation:  the  introduction  of  additional  equipment  into  a system 
to  reshape  its  root  locus  in  order  to  improve  system 
performance. 
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IV.  Control  System  Model  and  the  C Concept 

The  purpose  of  this  chapter  is  to  present  the  development  of  a 
discretely  controlled  system  model  of  the  ■YF-I6  incorporating  in  its 
design  consideration  the  handling  qualities  response  criteria  C*. 

First,  the  general  form  of  a discrete  servo  control  system  model  is 
presented.  This  is  followed  by  a discussion  of  the  C concept  and 
its  application  to  the  proposed  model. 

The  preceding  chapter  pointed  out  the  unstable  nature  of  Mach  .8 
flight  at  sea  level  for  the  YF-I6.  A requirement  for  some  type  of  con- 
trol system  for  this  novel  aircraft  in  this  flight  regime  was  repeatedly 
demonstrated.  At  this  junctiire,  several  options  are  available: 

a.  Synthesize  a continuous  control  law  using  classical 

, techniques  and  adapt  it  to  a digital  computer; 

b.  Synthesize  a control  law  using  continuous  optimal  control 
theoiy  and  again  adapt  it  for  a digital  computer;  or 

c.  Synthesize  a conti*ol  law  using  discrete  regulator/ 

• servo  theory  (Ref.  18) 

The  latter  approach  is  addressed  here,  due  to  its  direct  digital 
design  nature.  The  theoretical  approach  which  is  implemented  was 
first  developed  by  Sandell  (Ref.  12)  and  is  similar  to  the  more  recent 
proposeLLs  of  Stengel  (Ref.  15)  and  Lee  (Ref.  8).  Lee  points  out  the 
equivalence  of  his  approach  with  that  of  Stengel;  however,  Lee's 
application  is  to  a statically  stable  aircraft  for  a fixed  sample 
rate.  Here,  a discrete  optimal  control  formulation  with  a discrete 
quadratic  performance  index  will  be  tested  on  a statically  unstable 
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airframe  for  a variety  of  sample  rates  taking  servo  saturation  effects 
into  consideration. 

As  opposed  to  a regulator  problem  where  perturbed  states  are  re- 
turned to  zero  steady  state  values,  this  controller  scheme  will  produce 
a sequence  of  controls  to  force  the  trajectory  of  the  aircraft  to  follow 
some  input  reference  signal  introduced  by  the  pilot.  Pilot  dynamics  are 
not  considered  beyond  the  intimation  that  maneuvering  commands  originate 
with  him  with  no  time  delays. 

Control  System  Model 

Recall  that  the  continuous  control  system  model  as  depicted  in 
the  state  space  representation  of  Figure  17,  is  based  upon  a continuous, 
linear,  differential  system  of  equations.  The  system  is  considered 
deterministic  with  no  imcertainty  or  random  inputs  present  vdiile  u, 

X and  y,  represent  the  control,  state,  and  output  vectors,  respectively. 


It  should  be  noted  that  the  obseirver  matrix  with  its  associated  output 
vector  y could  be  incorporated  into  the  controller.  However,  for  the 
purpose  of  clarity,  it  is  desirable  to  separate  the  two.  The  above 
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figtire,  closely  represents  a regxilator  controller  in  that  no  external 
commands  enter  into  the  control  scheme.  If  €ui.  external  command  signal 
is  introduced  into  the  system,  the  regulator  problem  is  biased  to  fol- 
low that  signal.  The  problem  would  then  enter  the  realm  of  a servo  or 
tracking  problem.  Such  a situation  is  depicted  in  Figiire  18. 


Since  a commitment  to  a discrete  approach  for  digital  computer 
implementation  has  been  made,  a digital  equivalent  of  Figure  18  is 
necessary. 

The  proposed  equivalent  control  system  representation  appears  in 
Figure  19. 


The  symbolic  T in  Figure  19  represents  the  sample  rate  or  how 
often  samples  of  the  continuous  system  are  taken.  The  dashed  line 
indicates  that  all  switches  close  in  unison.  The  continuous  air- 
craft dynamic  model  is  represented  by  the  familiar  linear  equation 
set: 


Ax  + Bu 

(77) 

Cx 

(78) 

vrtiere  A,  B,  and  C are  the  dynamics,  control  and  observer  matrix, 
respectively  (Ref.  3).  Other  than  determining  what  A,  B,  and  C 
are,  it  is  necessary  to  define  the  system  states  (x),  the  system 
output  (y),  and  the  nature  of  the  input  command 

C Concept 

In  1966,  Tobie,  Elliott,  and  Malcom  introduced  the  C*  concept  as 
a new  longitudinal  performance  criterion  as  an  outgrowth  of  the  Cornell 
Aeronautical  Laboratory  "Thtunbprint"  (Fig.  26),  and  in  response  to  an 
effort  of  determining  vdiat  levels  and  types  of  handling  qualities  are 
required  by  pilots  (Ref.  17).  The  C criterion  is  an  example  of  speci- 
fying short  period  handling  qiialities  in  terms  of  aircraft  parameters 
familiar  to  a pilot  (6  aind  N^)  while  still  including  the  traditional 
short  period  frequency  (h.)  and  damping  reqirirements  ( ^ ).  As  these 
authors  point  out: 

"It  appears  likely  that  the  pilot  responds  to  the 
motions  that  naturally  tend  to  dominate  the  air- 
craft's characteristic  response.  For  example,  at 
low  velocity  where  cues  are  weak,  pitch  cues 
Would  be  most  importaint.  At  high  velocity  where 
very  slight  pitching  may  accompany  sizable  accelera- 
tion changes  Ng  cues  would  predominate.  It  is 
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postulated  that  the  pilot  then  responds  to  a 
blend  of  pitch  rate  and  normal  acceleration, 
with  the  blend  rates  varying  in  accordance 
with  natural  variations  in  aircraft  response. 
The  blend  of  response  parameters  has  been 
named  C*'." 


C is  defined  by  the  following  expression: 

C*  = K,  N„  + K.  9 + K~  5 

Z Z Q Q 

where  the  lonits/values  of  Table  IV  apply. 
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Table  IV 

C*  Terms  Defined 

Term 

Units 

Value 

C* 

g's 

variable 

g's 

variable 

& 

Rad/sec 

variable 

•• 

0 

Rad/sec^ 

variable 

Kz 

2 

g*  3-sec 

1.0  (arbitrarily  chosen  as  xmity) 

* 

g*  s-sec 

*co  (vdiere  V^q  = ADO  ft/sec) 

32.2 

"5 

g' s-sec 

distance  to  c.g.  from  pilot  _ 
gravity 

= .36 

32.2 

* 

Vgo!  cross  over  velocity,  at  which 
the  normal  acceleration  and 
pitch  rate  give  equal  cues 
to  the  pilot. 

50 


The  acceptable  range  of  the  C transient  response  to  a step  input  falls 
within  the  appropriate  bounds  of  the  C time  history  envelopes  of 
Figure  20.  The  four  regions  indicated  have  the  following  significance: 

I - Optijnal  response  (aerial  combat,  groimd  attack, 

and  penetration); 

II  - Not  as  critical  response  area  (air  refueling,  cruise) 

III  - Category  for  conditions  not  covered  by  I,  II  (loiter) 

IV  - Power  approach  (landing). 

These  regions  have  a direct  correlation  with  the  category  system  of 
MIL-F-8785B  (Ref.  9). 

Vlith  this  background  in  the  evaluation  and  meaning  of  the  C*  con- 
cept as  a well  defined  criterion  on  the  handling  qualities  of  an  air- 
craft, the  system  model  (Fig.  19)  can  be  put  into  perspective. 

In  this  study,  the  C envelope  will  be  used  to  measure  the  response 
of  the  system  with  satisfactory  response  defined  as  falling  within  the 

innermost  C*  boundary  limit.  The  output  of  the  system  (y(t))  is  re- 

. ^ ^ ^ 
placed  by  the  term  Ca,ct(t)  the  actual  system  C response. 

In  this  manner,  the  C*  response  becomes  a fiinction  in  the  perfomance 
index  of  the  system.  Additionally,  the  input  to  the  system 
is  replaced  by  which  represents  the  desired  C response  re- 

quested by  the  pilot.  In  a physical  sense,  the  YF-16  model  is  asked 
to  track  a pilot  step  input  command,  here  chosen  as  unity  since  it 
lies  in  the  center  of  the  envelope  boxmds.  In  the  YF-16,  using  a side 
motait  control  stick,  the  pilot  inputs  a desired  steady  state  value  of 

the  C parameter  to  the  controller  which  is  then  required  to  have  a 

* 

short  period  response  which  tracks  the  system  C response  to  the 
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commanded  C*  value.  It  is  apparent  that  the  observer  matrix  (C)  must 
transform  the  system  states  (x)  into  a blend  of  system  states  now 
redefined  by  equation  (79)  as  . 

The  next  chapter  will  determine  the  nature  of  this  transformation 
matrix  based  upon  a definition  of  the  system  states  in  equation  (77). 


E 
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^ Continuous  to  Discrete  System  Transformation 

This  chapter  presents  the  transformation  from  a continuous  sytem 
representation  to  a discrete  system  representation  necessary  to  continue 
develojment  toward  the  discrete  model  presented  in  Figrore  19. 

The  chapter  begins  with  the  continuous  system  state  model.  Here, 
the  nature  of  the  observer  matrix,  which  transforms  the  system  states 
into  a C representation,  is  developed.  The  chapter  concludes  with  a 
brief  discussion  of  the  discretization  process  used  to  represent  the 

i 

continuous  system  in  discrete  time. 

Continuous  System  State  Model 

Recall  from  Chapter  IV  that  an  open  loop  linear  continuous  system 
can  be  defined  by  the  first  order  state  variable  equation  set; 

— ~ x(t)  = 1^(T)  + B u(t)  (77) 

y(t)  = Cx(t)  (78) 

It  is  necessary,  therefore,  to  transform  the  three  original  equations 
of  motion  (equations  (l),  (2),  and  (3))  into  such  a representation  for 
their  implementation  in  the  controlled  system  model  developed  in  the 
previous  chapter  (Fig.  19).  This  is  done  as  detailed  in  Appendix  C 
using  the  short  period  approximation  with  the  resiilting  state  variable 


'=} 


•i 

• • 
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-2.603975 

15.058542 

0 


1.0 

-2.682339 

0 


-.260965 


-20.0 


u 

” 0 ” 
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& 

+ 

0 

•J 

A 

20.0 

(182) 


The  resulting  dynaniics  matrix  (A)  is  (3  x 3)  while  the  control 
matrix  (B)  is  seen  to  be  (3  x l) . A check  of  the  controllability  of 
the  system,  determined  by  assuring  that  j B]ABJ  A^  1?^  0,  shows  the 
system  is  completely  controllable. 

Referring  now  to  the  output  equation; 


y(t)  = C x(t) 


(78) 


We  can  replace  the  output  vector  y(t)  by  Cg^Q^(t)  as  discussed  in  the 
previous  chapter.  The  resulting  output  equation  is: 


Cact(^)  = ^ 


(80) 


As  shown  earlier,  C*  is  a linear  blend  of  3^  0 expressed  as; 


(79) 


C’'  = K^Na  + + K-0 

or  equivalently  as: 


pi  ''z  'S)] 


0 


(81) 


It  is  evident  that  a mismatch  exists  between  the  states  required 
to  calciilate  C , indicated  in  equation  (8l)  and  the  actual  system 
states,  developed  in  Appendix  C,  and  shown  in  equation  (18^.  Some 
transformation  matrix  D must  be  found  to  transform  system  states  to 
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the  required  states  needed  to  calculate  C . This  situation  is  ex- 
pressed in  equation  (82)  where  some  conformal,  (3  x 3),  D matrix  is 
shown  accomplishing  the  transformation. 


0 

Dll 

^12 

“13' 
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Nz 

= 

D2I 

D22 

D23 
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d 

•• 

© 

D31 

D32 

D33 

Cd) 

1 

(82) 

Three  linear  combinations  of  , & , and  must  be  fovind  which 

• •• 

are  equal  to  6,  N^,  and  0,  respectively. 

The  first  expression  is  the  trivial  relationship: 

d = + D]j2  0 + Djj 

0 = 0 + 1.0  © + 0 

For  the  next  equation, 

+ Doo  0 Doo  ^ 


(83) 

(84) 


= °21®^  °22  ® ^ ^23*^^ 


(85) 


recall  that. 


From  Appendix  C,  equation  (l74)  shows  ei  as: 

oC  = + 1.0  6 + 5. 

K,  Jf, 


(86) 


(174) 


This  equation  can  be  substituted  Into  equation  (86 ) leading  to.  the 
following  development: 
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^23  ~ 


i. 


(91) 


Similarly,  Appendix  C,  equation  6.78  ) shows  the  relationship  of  Q to 
oC  > 6 , and  as: 


^ (178) 


which  completes  the  development  of  the  D matrix  by  defining  D^l*  ®32> 


and  as:  j 

X 


c,  A ♦ c 
, ''  k^[  i y 


and,  _1_ 


Jti  . dzr  * 

k,  •<  K \ 


respectively. 
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Having  determined  D,  equation  (82)  can  be  vrritten  as: 


1.0  1 


N - 


1““— 

I k,  *w 


Substituting  this  expression  into  equation  (8l),  the  resulting 


expression  for  C is: 


This  expression  can  be  simplified  to: 


I * I 

. I 


Siibstituting  into  this  expression  the  appropriate  values  from  Table  III, 

•h 

Table  IV,  and  the  value  of  = 27.75155  based  on  ,8  Mach  at  sea 

W 

level,  re  stilts  in  the  following  expression: 


77.685  11.423  -9 


.921  . e 

L \ 
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i Equation  (95)  fixes  the  value  of  the  C obseirver  matrix  indicated  in 

Figure  19.  It  is  this  matrix  which  will  transform  the  system  states 
I into  the  system  scalar  output  C . Additionally,  the  system  is  found 

I to  be  completely  observable  since  the  determinant  of  |C  lA  C '(A  J C J 

is  nonsingular. 

The  continuous  system  described  by  equations  (l82)  and  (95)  is 
then  completely  controllable  and  observable  which  provide  the  sufficient 
! conditions  for  the  regulation  of  the  system  output. 

i 

Discrete  State  Model 

Because  of  the  need  to  implement  the  control  system  in  a digital 
flight  computer,  equation  set  (77)  and  (78)  must  be  transformed  into 
an  equivalent  discrete  form.  When  the  state  model  of  a continuous 
system  is  given  by  the  equation  set  (77)  and  (78),  the  discrete  state 


model  for  the  same  system  with  a piecewise  constant 

a difference  equation  set  of  the  form: 

input  is  given  by 

3E(K+1)T  = x(KT)  + u(KT) 

(96) 

y(KT)  = x(KT) 

(97) 

vrtiere, 

T at 

Ad  = ® 

(98) 

AT 

= oJ®  ^ (99) 


and  T is  the  sample  rate  (Ref.  2). 

Nxanerous  methods  are  available  for  the  calculation  of  the-  state 
transition  matrix  A^j.  For  example,  the  "Solution  by  Fimctions  of  a 
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Matrix"  method  (Ref.  2)  could  be  used  vri.th  the  values  of  the  and  B(j 
matrices  determined  by  hand  as  a general  function  of  the  variable  T. 
Additionally,  digital  computer  software  packages  are  available  to 
compute  the  values  of  and  . This  is  the  method  chosen  in  this 
study  since  the  matrices  may  be  determined  rapidly  for  any  specified 
sample  rate  in  a call  to  a subroutine  nested  in  a more  complex  main 
program. 

In  particular,  and  B^  are  digitally  calculated  in  this  re- 
search by  use  of  the  digital  subroutine  DSCRT  discussed  in  Ref.  6. 

In  the  calculation,  B^  is  determined  by  e 'aliiating  the  series: 


5<i  = H ^ ^ * 


vrtiile, 

A^  = I + A . Bjj 

trtiere  ^ and  NT  are  specified  in  the  call  to  the  subroutine 
(see  Appendix  E). 

Examples  of  k^  and  B^  for  vjirious  values  of  T are  listed  in 
Table  V. 


(lOO) 


(101) 
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Table  V 


Ajj  and  Matid-ces  as  a Fimction  of  Sample  Rate 


T = 1/50  sec 

Bd 

• 

.9521122 

.2859504 

0 

.0189892 

.9506241 

0 

-.0122758 

-.7653255 

.6703200 

-.0020253 

-.1647557 

.3296800 

T = 1/70  sec 

.9650 

.2073 

0 

.01376 

.9639 

0 

-.997504 

-.5814 

.7515 

-.0009032 

-.08757 

.2485 

T = 1/100  sec 

.9750 

.1467 

0 

.009742 

.9743 

0 

-.004528  ^ 
-.4265 
.8187 

-.0003916 

-.04427 

.1813 

r 

I 

f 

! VI.  ZOH,  Optimal  Controller  and  Simulation  PeveloFment 

This  chapter  presents  the  development  of  the  sequence  of  optimal 
scalar  controls  u(KT)  which  will  drive  the  system  shown  in  Figure  19 

i * 

' to  the  commanded  steady  state  value  of  C input  by  the  pilot.  As 

Figure  19  indicates,  the  controller  accepts  sampled  information  about 
the  system  and  uses  this  information  to  construct  control  inputs  to 

the  system.  The  determination  of  the  control  law  will  be  accomplished  I 

; I 

by  minimizing  an  appropriate  discrete  penalty  function  with  zero  steady  i 

state  error.  The  resulting  controls  are  step-like  in  nature  as  shown 
in  Figure  21. 


Fig.  21.  Nat\ire  of  Optimal  Controls. 

From  Figure  21,  it  is  apparent  that  the  control  u(KT)  takes 
some  constant  value  during  the  time  period  between  samples.  That  is: 


u(KT)  = constant  V (KT)  ^ -f  < (K+1)T,  K = 0,1,2,3,...N 

(102) 
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Additionally,  the  sample  interval  (t)  remains  constcint  while  the  magni- 
tude of  the  control  is  allowed  to  vary.  This  series  of  piecewise  con- 
stant step  controls  corresponds  to  the  shift-like,  control  stick  input 
commands  introduced  by  a pilot.  Altho\agh  a step  command  may  not  re- 
present the  most  common  pilot  input,  the  step  response  implicitly 

I 

describes  the  response  to  other  inputs  that  a pilot  uses.  Additionally,  j 

Figure  21  is  indicative  of  the  output  of  a Zero  Order  Hold  device.  | 

This  chapter  will  first  discuss  the  general  concepts  of  a hold 
device  and  in  particular,  the  Zero  Order  Hold.  Next,  a discussion  of 
the  discrete  penalty  function  selected  with  respect  to  which  the  per- 
formance of  the  system  is  optimized  is  presented.  The  recursive  con- 
trol algorithm  which  minimizes  this  cost  function  is  then  incorporated 
in  a computer  program  developed  for  this  investigation.  A discussion 
of  the  structiiring  of  this  simulation  program  is  presented.  The  chapter 
concludes  with  the  development  of  the  closed  loop  system  eigenvalues. 

The  determination  of  these  eigenvalues  is  incorporated  as  a part  of 
the  softwaire  and  used  to  track  the  migration  of  system  roots  as  a 
function  of  sample  rate. 

Zero  Order  Hold  Device 

The  type  of  hold  device  incorporated  in  a digital,  control  system 
cam  play  a decisive  role  in  determining  vdiether  the  system  is  stable 
especiaiUy  at  lower  sample  rates,  A hold  device  serves  to  convert  a 
discrete  time  sequence  of  numbers,  sepairated  in  time  by  T-second 

intervads,  into  a continuous  time  function  in  order  to  provide  a ! 

suitable  input  to  a continuous-time  component.  The  pamticulao*  hold  j 

device  chosen  represents  a D/A  converter  which  constructs  a piecewise  ^ 


k: 
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continuous  control  signal  from  a pialse  sequence  of  numbers  each  of 
which  represent  a discrete  control  signal.  For  a control  problem 
such  as  this,  it  is  more  appropriate  to  think  in  terms  of  extrapolating 
the  present  control  over  the  time  interval  between  samples  rather  than 
trying  to  reconstruct  some  signal  which  is  more  along  the  vein  of  a 
communications  problem. 

Between  sampling  instances,  the  hold  device  extrapolates  between 
the  most  recent  sample  and  the  next  to  follow. 

A power  series  expansion  of  a continuous  signal  u(t)  in  the 
interval  between  sampling  instant  KT  and  (K+l)Tcan  be  expressed  as; 

u(t)  = u(KT)  + d [ u(KT)  ] (t-KT)  + d^  [u(KT)1  (t-KX)^  + ... 

dT  1!  dT  2! 

VKTSt<^K'»0T  (103) 

Using  the  higher  order  derivatives  of  u(t),  for  the  purpose  of  more 
accurate  extrapolation,  can  meet  with  serious  diffic\alties  in  main- 
taining system  stability.  This  is  due  to  the  fact  that  the  higher 
the  order  of  the  derivative  to  be  approximated,  the  larger  the  number 
of  delay  pxolses  required  since  one  time  interval  must  pass  for  each 
discrete  control  needed  in  the  power  series.  The  accuracy  of  the 
estimate  of  a derivative  is  then  a function  of  the  number  of  time 
delays.  It  is  these  time  delays  which  have  an  adverse  effect  on  the 
stability  of  a feedback  control  system.  The  value  of  the  first  deriva- 
tive is  known  from  the  calcvLLus  to  be  approximated  by  the  simple 
difference  equation; 

- u(t)  ( 10k) 


1 

T 


u(KT)  - u(K-l)T 
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It  is  apparent  that  for  this  first  derivative  approximation,  both  the 
present  sample  value  and  the  immediate  past  sample  value  are  needed. 

For  higher  order  derivatives,  the  number  of  past  control  pulses  ,, 

li 

required  grows  larger  and  larger.  For  example,  the  second  derivative;  jj 


I 

I 

I 

I 

i 

I 


u(t)  - 


2^ 

rp2 


u(KT)  - 2u(K-l)T  + u(K-2)T 


(105) 


requires  three  consecutive  controls.  Equation  (103)  can  be  thotight 
of  as  a best-fit  m^^  order  polynomial  approximation  of  u(t)  which 
may  be  rewritten  as; 


Ult)  = u(KT  + T ) = % A.T  ^ A,  ( 106) 


where: 
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vrtiere  is  equated  to  u(KT).  This  makes  sense  since  it  is  only 

natural  to  require  the  output  signal  u(t)  to  have  the  value  of  the 
input  sequence  at  the  sample  time.  The  Zero  Order  Hold  is  therefore 
expressed  as; 

u(KT  + T ) = u(KT)  (109) 

As  Figure  21  confirms,  the  effect  of  the  ZOH  is  to  stretch  the 
input  pulse  into  a series  of  rectangular  waves  of  width  T.  Finally, 
it  should  be  noted  in  the  change  in  system  representation  from  Figure 
19  to  Figure  23  that  the  ZOH  effect  is  taken  into  account  in  the  pro- 
cess of  discretizing  equation  (77)  into  equation  (96). 

Penalty  Function 

The  cost  criterion  selected,  with  respect  to  vdiich  the  perfor- 
mance of  the  system  is  optimized,  is  of  the  form  of  the  quadratic 
functional; 

u '^(t)  R u(t)  ( dt  (no) 

J (Ref.  8) 

where  Q is  the  weighting  which  penalizes  the  trajectory  deviation  or 
difference  between  the  C commanded  and  the  C realized,  i.e.,  Cx(t), 
and  R is  the  weighting  on  the  control  vdiere  a penalty  is  exacted  for 
larf*  corrective  control  rates. 

It  Is  required  to  determine  a discrete  control  sequence  vdiich 

eifiia:**  • Kecret.e  oquivalent  of  equation  (lio)*  Using  the 

f 'be  first  derivative  of  the  control 


u(t),  the  cost  functional  (J)  of  eqviation  (llO)  can  be  replaced  by 
the  following  equivalent  discrete  expression: 


^d  =E  ( 

ko»  - = 

T 

Qd 

Cm  - C x(KT)  1 

k»ol 

L J 

u(K+l)T  - u(KT) 

T 

^d 

-• 

u(K+1)T-  u(KT)jj 

(in) 


or  equivalently  as: 


- f { [4.  - =:=t] " 


K = 0 

Ju(K+l)T  - u(KT)j  ^ J 


I u(K+l)T  - u(KT) 


] 


(1S5) 


As  pointed  out  in  Appendix  D,  the  control  which  minimizes  this  cost 
function  is  given  by: 

u(K+l)T  - u(KT)  = - cjct)  + 


Nd  [ x(K+l)T  - x(KT)]  {235) 


which,  when  evaluated  for  a few  values  of  K,  i.e,,  K = 0,  1,  2 . . 
can  be  expressed  by  the  recursive  equation: 


K-1  ^ 1 

u(KT)  = X;  ^d  ^Cm  - <ct)  '‘d  [ - x(0)J+  u 


(0) 


J=0 


(112) 


where  x(0)  and  u(0)  are  the  initial  states  and  control  and  where: 


_ _ _ -1  _ _ _ _ -1  _ -1 

Id  = ^^2.  - Ki  (A^  - I)  )(Cd  (Ad  - I)  Bd) 

d d 


(:?38) 
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(240) 


"d  = ('i,  + h R 

d 

The  matrices  K,  , Ko  are  foimd  from  the  positive  definite  solution 
■^d  d 

of  the  algebraic  matrix  Ricatti  equation  of  the  fom: 

- Wtu-"  Tfi  (219) 

(Ref.  6) 

or  from  Appendix  D for  this  particular  problem: 
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(113) 


where  the  resiilting  steady  state  solution  matrix 


used  in  the  gain  calculation  equation  expressed  by: 

-[r^r.  ./rt'i 


Pn 

d 

P12 

( 

P « 

12d 

22^ 

(Ref.  8) 

d 

d 


is 


(218) 


or  more  precisely  as: 
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[21] 

“ — 1 
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+ 

• [-] 

-lid  -12d 

h §d 
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- 
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(UA) 


where  Pn  , P , and  P„  for  this  application  are  matrices  of 
d -12(i 

dimension  (3  x 3),  (3  x 1),  and  (1  x 1),  respectively. 

Computer  Program  and  Simulation 

.»•  Eq\iation  (112)  is  suitable  for  recursive  evaluation  on  a digital 

computer,  along  with  the  terms  which  comprise  it  expressed  by  equations 
(219),  (238),  and  (240).  Appendix  E,  contaiins  the  digital  com- 
*'■>,  puter  program  which  follows  the  analytical  solution  to  the  servo 

problem  expressed  by  the  previous  set  of  equations  and  simulates  a 
system  eqvdvalent  to  Figxire  19.  A simplified  flow  diagram  of  the 
program  is  shown  in  Figure  22.  The  structure  of  the  simulated  closed- 
loop  system  appears  in  Figure  23  while  its  programming  logic  appears 
in  Figure  24.  The  linear,  tirae-invarient,  continuous-time  system 
given  by  equations  (77)  and  (78),  having  been  transformed  into  the 
equivalent  discrete-time  linear  equation  set; 


x(K-H)T  = x(KT)  + Bjj  u(KT) 

(96) 

C*  (KT)  = x(KT) 

(97) 

where  ~ ^ > with  a quadratic  cost  functional  used  to  determine 
u(KT),  is  simxilated  in  the  algorithm.  All  the  required  matrices  are 
sequentially  evaluated  for  a specific  sampling  interval  with  the 
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START 


Specify  Sample 
Rate  (TDEL) 


Set  Simulation 

Flag  — >-(Flag  other  than  1 

1 discussed  later) 


Specify  Simulation 
duration  in  seconds 
(LONG) 


Adjust  sample  rate  for 
computer  simulation 
use  (del) 


Load  in  continuous  system  A,  B,  and  C 
matrices.  (Built  internally  into 
program  for  this  one  flight  condition.) 


Specify  Q and  R 
weighting  matrices 


Calculate  ly  and  N^j 


Run  ZOH  Simulation 


List  Results 


CalcTilate  System  Roots 


Plot  Results  optional 

(See  Appendix  F) 


Fig.  22.  Simplified  Computer  Program  Flow  Diagram. 
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Simulation  Closed  Loop  System, 


start  Simulation 


Initialize  all 
Variables 


Print  Listing  and 
Write  on  Tape  Initial 
Values 


Increase  Time  and 
Iteration  Counters 


^^Time  N. 
for  controP 
Update?/” 


Calculate  New 
States  and  C* 
Output 


Print  Listing  and 
Write  on  Tape  New 
Values  of  states, 
C*,  and  control 


“^to  Stop^ 
Simulatior 


Increase  Time  and 
Iteration  coiaiters 


Calculate  new 
Control 


Calculate  new 
states  and  C'*' 
output 


Fig.  2L.  Simplified  Simulation  Flow  Design. 
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RLcatti  matrix  equation  steady  state  solution  used  to  calculate  the 
matrix  of  feedback  gains.  The  ailgorithm  was  programmed  in  Fortran  IV 
and  was  tested  on  a CDC  6600  digital  computer. 


Insert  (a)  of  Figure  23  is  a breakdown  of  the  discrete  plant 
dynamics.  Since  the  system  being  modeled  is  continuous,  some  way 
must  be  found  to  simulate  it  on  a digital  computer.  For  this  problem, 
the  continuous  system  was  discretized  at  a rate  of  1/500^^  second. 

This  rate  of  plant  discretization  is  kept  constant  no  matter  what  the 
sample  rate  specified  for  investigation  might  be.  In  other  words,  the 
YF-16  plant  is  allOTired  to  change  states  every  1/500^^  second.  Such  a 
fast  rate  allows  the  plant  to  retain  an  almost  continuous,  real-world, 
nature  while  still  satisfying  the  requirement  for  discretization.  It 
is  important  to  realize  that  the  sample  rate  specified  for  a particular 
program  run  is  the  discrete  rate  with  vdiich  this  equally  discretized 
plant  model  is  sampled. 

Initially,  the  plant  was  discretized  to  a rate  of  some  multiple 
of  the.  sample  rate.  For  example,  if  the  sample  rate  was  set  to 
T = .333  sec.,  the  plant  was  always  discretized  to  a rate  ten  times 
faster  or  .033  seconds.  However,  to  standardize  the  runs,  thus 
allowing  a comparison  of  results,  and  remove  this  factor  of  ten, 
which  might  influence  the  results  it  was  decided  to  standardize 
the  plant  discretization  rate.  With  the  "factor  of  ten"  scheme, 
however,  it  was  much  simpler  to  determine  vdien  it  was  time  to  bxiild 
a new  control.  For  example,  if  the  sample  rate  specified  was 
T = 1/40^^  second,  ten  iterations  of  the  plant  ( ~ ), 

this  time  running  at  l/400^^  second,  passed  before  a new  control 
was  calculated.  This  latest  control  was  then  implemented  and  retained 
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as  the  plant  control  for  the  next  ten  plant  iterations  (ZOH  effect) 
before  the  cycle  was  repeated.  It  is  evident  that  this  factor  of  ten, 
while  simplifying  the  iteration  and  update  scheme,  made  the  correlation 
of  resiilts  from  one  sample  rate  to  another  almost  impossible. 

The  standardization  scheme,  while  eliminating  this  problem, 
introduced  the  new  problem  of  how  to  determine  when  to  update  the 
control.  The  frequency  of  control  updates  now  varied  with  the  re- 
lationship  of  the  sample  rate  to  the  l/500^  second  plant  discretization 
rate.  If  a sample  rate  of  T = l/25  second  was  used,  twenty  plant 
iterations  passed  (20/5C)0  = l/25  second)  before  control  updates. 

This  was  simple  enough  for  even  multiples  of  the  plant  discreti- 
zation rate  but  more  complicated  for  non^nultiple  cases.  The  problem 
could  be  eliminated  by  specifying  the  use  of  only  convenient  multiple 
sample  rates  (i.e.,  T = V^O,  1/50,  l/lOO,  etc,).  However,  to 

rotain  program  flexibility  in  the  use  of  any  sample  rate,  this  limi- 
tation was  not  implemented.  Instead,  the  flexibility  was  retained 
but  at  .the  cost  of  some  slight  error.  The  Fortran  function  IFU  is 
xised  in  conjunction  with  the  modulo  arithmatic  function  MOD.  This  is 
best  explained  by  use  of  an  example.  If  the  sample  rate  is  specified 
to  be  T = 1/30^^  second,  ^rtiich  in  the  program  is  referred  to  as 
TDEL  = .0333...,  the  program  used  the  Fortran  expressions: 

MODI  *=  IFU  (500.0  * TDEL  + .5) 

DEL  *=  MODl/500.0 

to  come  up  with  the  adjusted  sample  rate  of  DEL  = ,034  seconds. 

The  value  of  DEL,  in  this  case  .034  seconds,  is  then  used  through- 
out the  remainder  of  the  program  as  the  adjusted  sample  rate  to  simplify 
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the  calculations.  The  amount  of  error  introduced,  in  this  case,  the 
difference  betvreen  .0333...  and  .034  seconds,  is  negligible.  A com- 
parison of  the  sample  rate  specified,  to  the  actual  rate  selected  by 
the  algorithm  for  the  simulation  appears  in  Table  VI.  It  is  evident 
that  for  specified  sample  rates  , vrtiere  500/r  equals  an  integer 
exactness  between  the  sample  rate  specified  and  the  sample  rate  used 
is  maintained. 


Specified 

Table  VI 

vs  Actual  Sample  Rates 

Specified  Sample 

Actual  Sample 

Rate  (del) 

Rate  Used  (TDEL) 

1/100.0 

1/100.0 

1/90.0 

1/83.3 

1/80.0 

1/83.3 

1/70.0 

1/71.4 

1/60.0 

1/62.5 

1/50.0 

1/50.0 

1/40.0 

1/38.4 

1/30.0 

1/29.4 

1/20.0 

1/20.0 

1/10.0 

1/10.0 

Additionally,  the  value  of  MODI  is  \ised  to  determine  when  a control 
update  should  occiir  using  the  Fortran  modulo  (MOD)  arithmetic 
statement : 

(MOD  (KK,  MODI) J 2,  1,  2 

Here,  KK  is  an  integer  counter  vrtiich  keeps  track  of  the  number  of 
system  iterations.  If  the  value  of  KK  is  divisible  by  MODI  such  that 
the  remainder  is  zero,  the  program  will  jump  to  statement  1 vdiere  an 
update  of  the  control  occurs.  For  any  other  value  of  remainder,  the 
program  will  continue  using  its  present  value  of  control.  Such  a 
scheme  was  used  successfully  in  this  study. 

Insert  (b)  of  Figure  23  is  a detaiiled  depiction  of  the  mechani- 
zation of  the  discrete  controller  expressed  by  equation  (96).  The 
optimal  scalar  gain  and  the  (1x3)  optima]  gain  matrix  are 
calculated  based  on  the  modified  sample  rate,  DEL,  just  discussed. 
Since  the  control  is  calculated  every  DEL  seconds,  it  is  necessary 
that  the  gains  be  based  on  this  same  rate.  The  values  of  L^  and 
are  determined  by  sequentiail  calls  to  subroutines  as  shown  in  Figure 
25.  Additionally,  since  the  input  to  the  simulation  is  a step  of 
unity  magnitude,  a sample  taken  at  any  time  will  also  be  unity. 

That  is: 

Cqqju  = ^com  ~ 1«0.  (^5) 

discrete  continuous 

Finally,  with  this  unity  step  input  driving  the  controller,  and 
assuming  the  initial  control  u(0)  aind  initial  states  x(0)  are  taken 

as  zero,  it  is  apparent  that  the  first  control  calculated  will  have 
the  value  of  L^.  This  is  the  approach  taken  in  the  simulation. 
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From  Main  Program 


(CALL  DSCRT) 


Discretize  continuous  A & B 
matrices  to  DEL  sample  rate 


(CALL  RIC) 

Detennine  steady  state  solution 
to  Ricatti  matrix  equation 


(CALL  GAINK) 


Determine  K-.  (1x3) 

d 

and  Ko  (l  x l)  values 
^d 


(CALL  NDLD) 


Calculate  values  of 
(1  x 3)  and  ly  (1  X 1) 


GO  TO  SIMULATION 


Fig.  25.  Simplified  Flow  Diagram  of  I^  and  Calculation. 

A key  element  in  the  determination  of  the  control  u(KT)  is  the 
amount  of  computation  time  required  for  its  calculation.  With  the 
plant  running  at  a 1/500^*^  second  state-change  rate,  the  computation 
time  required  for  both  state-output  and  control  calculations  must  be 
much  less  than  this  rate.  This  is  especially  true  for  the  control 


computation  since  it  must  be  available  on  the  next  simulation  iteration 
which  occurs  l/5CXD^^  second  later.  Such  a state  of  aifairs  is  true  not 
only  for  this  simulation  but  becomes  even  more  critical  when  the  com- 
putations are  done  on  the  type  of  onboard  computer  envisioned  as  the 
actual  aircraft  controller.  Such  computers  are  of  much  more  limited 
capability  in  terms  of  wordlength  and  speed  of  computation  than  the 
ground  based  GDC  computer  used  here.  Algorithm  computational  time  can 
be  kept  at  a minimum  by  efficient  programming,  which  limits  the  number 
of  computer  commands,  and  by  proper  sequencing  of  these  computer  com- 
mands. The  most  time  consuming  between  iteration  computation  in  this 
simulation  is  taken  by  the  sequence  of  statements  which  calcxilate  the 
new  control.  The  calculation  process  has  been  reduced  to  five  sequential 
statements  by  retaining  information  from  previous  control  computations 
and,  for  this  reason,  is  assumed  to  occur  in  much,  much  less  than 
1/500^^  second. 

This  assmption  is  reasonable  in  view  of  the  amo\mts  of  time  small 
computers  require  for  computation  piirposes.  Working  on  the  micro-second 
(10“^  second)  level,  no  time  problems  are  envisioned  in  the  various  data 
manipulations  (shift  commands,  access  memory  location  commands,  comple- 
ment commands,  etc.)  pxirsuant  to  the  calculation  of  u(KT),  By  the 
elimination  of  the  use  of  intermediate  answers,  as  presently  programmed, 
memory  access  time  could  be  minimized  even  further  lending  greater 
credence  to  this  assumption. 

This  assiimption  of  negligible  lapse  in  time  for  computation  pur- 
poses is  of  additional  importance  vrtien  equation  (II2)  is  examined 

i 

closely.  From  this  eqviation,  it  is  evident  that  the  calculation  of  ’ 

a(KT)  requires  knowledge  of  x(KT),  The  x{KT)  used  is  the  most  ciirrent  ! 
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I set  of  values  available.  The  resvilting  u(KT)  is  actually  u(KT'*’)  owing 

to  a finite  amount  of  computational  time  delay.  This  uCKT"^)  is  then 
used  as  the  control  u(KT)  in  subsequent  state  evaluations  using 
equation  (96  ) . 

; Closed  Loop  System  Eigenvalues 

An  additional  feature  of  the  program  is  the  inclusion  of  a sub- 

^ \ 

routine  (ROOT)  which  calcxilates  the  characteristic  roots  of  the  closed- 
loop  system.  The  location  of  roots  inside  the  unit  circle  and  their 

i 

^ migration  with  changes  in  the  sample  rate  is  of  interest.  Subroutine 

7 

I ROOT  determines  these  roots  by  solving  for  the  eigenvalues  of  the 

i dynamics  matrix  of  an  augmented  state  space  system  representation. 


This  dynamics  matrix  is  determined 

as  follows.  Let 

K-1 

”1 

W(KT) 

S ^d 
i=o 

Ccom(KT)  - cIct(KT)J 

(116) 

then. 

1 

b\it. 

W(K+1)T  = 

W(KT)  + 

Okt)  - 4t(KT)j 

(117) 

- 

so. 

►a 

II 

Cd  x(KT) 

(97) 

or. 

W(K+1)T  = 

W(KT)  + ly 

(^com^KT)  - Cd  x(KT)) 

(118) 

W(K+1) 

^?(KT)  + 

Cgom(^T)  ~ ^ ^d 

(119) 
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This  last  eqiiation  characterizes  the  dynamics  of  the  closed-loop  system 
of  Figure  23.  The  eigenvalues  of  the  (2x2)  dynamics  matrix  of  this 
equation,  as  a function  of  the  sample  rate,  define  the  stability  of  the 
system  in  the  Z-plane,  Subroutine  ROOT  calcxolates  these  eigenvalues. 

It  is  desirable  that  the  resulting  short  period  roots  not  only  lie  with- 
in the  \aiit  circle,  insuring  system  stability,  but  aiso  that  they  corre- 
spondto  a naturad  frequency  ( and  damping  ratio  ( ) that  are 

preferred  by  pilots.  Preferred  values  of  verses  ^ are  given 
graphical  significance  in  what  is  called  the  Cornell  Aeronautical  Lab 
"Thumbprint"  (Fig.  26). 

This  particular  vs  ^ relationship  is  based  on  comprehensive 
data  from  severail  variable  stability  airplanes,  primarily  the  F-94. 

The  various  values  and  their  associated  ^ values  can  be  mapped 
into  the  S-plane  and  results  in  the  "kidney-shaped"  representation 
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of  Figxire  27.  In  the  Z-plane,  Figvire  26  corresponds  to  a similarly- 
shaped  region,  albeit  much  smaller,  shown  as  I in  Figure  28.  The 
exact  size  of  the  "thumbprint"  in  the  Z-plane  is  determined  by  the 
sample  rate  selected.  As  the  sample  rate  increases,  approaching  con- 
tinuous sampling,  i.e.,  T =1'"'  = 0,  the  Z-plane  contour  shrinks 

and  converges  to  the  single  point  at  u = + 1.  If  the  sample  rate  is 
reduced,  the  contour  expands  ^rtiile  at  the  same  time  moving  away  from 
the  u = +1  point  as  region  II  of  Figure  28  depicts.  If  the  center  of 
the  thumbprint  (Fig,  26)  is  chosen  as  the  most  desirable  operating 
point,  corresponding  to  $ = .707  and  the  following 

conjugate  set  of  S-plane  root  locations  resvilt; 

= -3*89  ± 3.389 j (liJ3) 

Using  the  direct  Z-transform,  Z = e , this  corresponds  to  pole 
locations  in  the  Z-plane  of; 

2 = ± 3.89j)T  (124) 


I 
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.Vj  Z = u + vj 

+1 


Fig.  28.  Thvunbprint  in  Z-plane. 

>rhioh  for  T = 1/50  becomes: 

-3.89  +3.89,1 

Z = e 50  * e 50 

= .925149  / + 4.457“ 

(125) 

= .925149  (cos  4.457  ± j sin  4.457 
Z = .92235  ± .0719j 

The  Z-plane  roots  resulting  from  subroutine  ROOT,  which  are  fixed  by 
the  optimal  solution  to  the  problem,  will  be  looked  at  with  respect 
to  this  most  desirable  root  location.  It  is  evident  from  equation 
(124)  that  the  particular  location  of  the  roots  will  vary  with  the 
sample  rate  selected. 


In  summary,  this  chapter,  which  begin  with  a brief  discussion  of 
hold  devices,  showed  the  development  of  a discrete  cost  criterion 
(penalty)  function  and  its  minimizing  solution.  A computer  program 
implementing  this  optimal  solution  and  some  of  the  problems  associated 
with  its  development  were  then  discussed.  The  chapter  concluded  with 
a discussion  of  the  determination  of  the  system  closed  loop  eigenvalues. 
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VH.  Slmxilation  Re  stilts 


This  chapter  presents  the  resxilts  of  the  digital  simulation.  The 
digital  program,  listed  in  Appendix  E,  was  exercised  for  various  sample 
rates  and  weighting  factors  (R  & Q)  with  the  following  prioritized 
questions  in  mind: 

a.  Can  the  system  be  controller? 

b.  Is  the  controlled  response  within  the  C* 
envelope  bounds? 

c.  Is  the  solution  a feasible  one? 

An  affirmative  answer  to  question  (a)  was  achieved  and  the  system 
was  considered  controlled  if  the  system  successfully  tracked  the  1-G 
climb  command  input  to  the  system.  This  was  evident  idien  the  steady 
state  error  between  the  C*  input  command  and  the  actual  system  C* 
response  was  zero.  Additionally,  in  such  a steady  state  condition, 
the  system  states  or  , 0 , and  6.  would  also  achieve  steady  state 

values.  A divergent,  undamped  response  as  Figures  11a  and  12a  depict, 
would  be  unacceptable  and  indicative  of  a total  lack  of  control. 

Having  achieved  a controlled  system,  question  (b)  coild  be  an- 
swered  by  comparing  the  system  C response  to  the  established  C 
envelope  bounds  of  Figure  20.  To  do  this,  a rescaled  plastic  overlay 
of  Figrore  20  was  superimposed  over  the  C response  for  a particular 
program  run.  In  this  way,  cases  \^ch  fell  outside  the  bounds  could 
be  identified. 

The  final  question  of  achieving  a feasible  solution,  cotild  be 
answered  by  observing  the  gains  (ly,  N^j)  determined  as  necessary  by 
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the  program,  and  the  control  deflections  they  produced.  A set  of  gains 
which  control  the  system  and  cause  the  response  to  fall  within  the  re- 
quired envelope  but  which  are  infeasible  in  terms  of  the  amount  of  con- 
trol sxtrface  deflection  they  require  for  implementation  or  the  rate  of 
control  application  (satioration  effect),  would  also  prove  unacceptable. 

The  initial  ZOH  simulation  runs  were  conducted  with  arbitrary 
values  for  both  the  Q and  R weightings.  After  considerable  trial  and 
error,  it  became  evident  that  a Q weighting  of  unity  produced  a consis- 
tent maxiraum  response  overshoot  of  the  target  value  (l-G)  of  approxi- 
mately five  percent.  It  was  convenient  to  think  of  the  Q weighting 
(trajectory  error  penalty)  as  determining  the  amount  of  overshoot  while 
the  R weighting  (energy  expenditure  penalty)  determined  the  amount  of 
control  deflection  employed. 

Subsequent  simulation  runs  were  therefore  conducted  keeping  the  Q 
weighting  equal  to  unity  while  the  sample  rate,  and  R weighting 
(affecting  the  ratio  of  Q to  R)  were  allowed  to  vary. 

Typical  Results 

ZOa  simulation  runs  were  reinitiated  with  both  the  trajectory  error 
weighting  and  control  penalty  weighting  set  to  unity.  These  values  were 
selected  to  obtain  an  initial  feel  for  the  system  response  and  some  idea 
of  the  nature  of  the  controls  produced.  The  program  output  for  a typical 
simulation  run  is  contedned  in  Appendix  F,  When  this  output  information 
is  made  available  to  the  plotting  program  listed  in  Appendix  G,  the  out- 
put products  of  Figures  29,  30,  and  31  result.  These  particular  plots 
are  for  a sample  rate  of  T = l/50^^  second  with  Q = 1,0  and  R = 1,0. 
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Figure  29  presents  a two  second  time  idstory  of  the  system  states 
6 , and  6,  . From  this  figure,  it  is  appareiit  that  after 

approximately  .80  seconds,  the  states  have  achieved  a steady  state 
condition.  The  initial  horizontal  stabilizer  deflection  is  negative 
which  causes  the  nose  of  the  aircraft  to  pitch  up  in  almost  ininediate 
response  to  the  climb  command  initiated  at  zero  seconds.  The  controller, 
sensing  the  reducing  magnitude  of  the  error  between  C actual  and  C 
command,  gradually  reverses  the  stabilizer  deflection  direction  thus 
lowering  both  the  pitch  rate  and  the  angle  of  attack  to  their  eventual 
steady  state  values.  In  the  steady  state,  a slightly  positive  deflec- 
tion of  the  stabilizer  is  required  to  maintain  the  climb. 

The  next  figure  in  the  series.  Figure  30,  shows  both  the  actual  C 
output  of  the  system  and  its  relation,  as  a f motion  of  time,  to  the 
commanded  response. 

After  some  initial  hesitation,  due  to  the  delay  before  a sample  of 
the  system  states  occurs,  the  controller  responds  rapidly  to  drive  the 
system  to  the  commanded  response  with  an  overshoot  of  5.1  percent  before 
reaching  a steady  state  lock-on  to  the  commanded  response  in  .556  seconds. 

The  final  figvire  in  the  series.  Figure  31,  shows  a time  history  of 
the  controls  which  produce  the  state  and  output  responses  just  discussed. 
Initially,  the  controls  remain  at  zero  until  enough  time  (sample  interval) 
has  passed  for  a sample  to  be  taken  and  a control  calculated.  In  this 
case,  this  occurs  after  .02  seconds  with  the  initial  control  equaling 
the  value  of  the  gain  The  control  surface  is  deflected  negatively 

a maximm  of  .01911  radians,  then  positively  a maximum  of  .005695 
radians.  It  finally  stabilizes  to  a positive  value  of  .001571  radians 
after  ,80  seconds  vrtiich  corresponds  to  the  same  elapsed  time  needed  for 
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the  states  to  stabilize.  The  Zero  Order  Hold  effect  of  maintaining  the 
controls  constant  between  samples  is  clearly  evidenced  by  those  portions 
of  the  plot  which  remain  constant  over  time  before  jtanping  to  a new 
constant  value. 

Modification  of  Control  Scheme 

In  an  effort  to  reduce  the  amoimt  of  time  it  took  for  the  actual  C 
response  to  achieve  and  maintain  C'**’  commanded,  with  zero  steady  state 
error,  shown  in  this  example  (Fig.  30)  to  take  .556  seconds,  a modifi- 
cation to  the  hold  technique  used  in  the  simulation  was  undertaken. 

From  observing  nvunerous  control  time  history  plots,  of  vdiich  Figure  31 
is  a typical  example,  it  was  felt  that  time  could  be  conserved  in  the 
control  scheme.  Conserving  time,  in  this  case,  could  reduce  the  time 
required  to  reach  steady  state  to  something  less  than  .556  seconds. 

After  the  calcialation  and  implementation  of  the  first  non-zero  control, 
for  example,  the  controls  remain  negative  and  constant  until  the  next 
update.  At  the  next  update,  it  is  apparent  from  Figure  31  that  a yet 
more  negative  control  results  from  the  update  computation.  It  wotild 
appear  advantageous  then  to  use  some  of  the  time  the  controls  are  held 
constant  to  move  toward  this  new  control.  This  would  appear  especially 
helpftil  in  saving  time  on  the  positive,  upward  sloped,  portion  of  the 
plot  (Fig.  31).  Here,  the  control  deflection  has  reversed  and  begins 
to  move  positive  until  reaching  its  peak  of  .005695  radians  at  .202 
seconds  elapsed  time.  The  scheme  of  the  modified  control  would  appear 
as  in  Figure  32. 

Fig\ire  32  is,  in  effect,  a First  Order  Hold  (FOH)  representation 
of  the  control,  vdiere  the  control  changes  at  a constant  rate  between 


91 


Fig.  32.  Updated  Control  Plus  Sloped  Projection  Scheme. 


discrete  update  intervals  (KT).  On  any  particular  iteration,  the  con- 


trol can  be  expressed  as: 


u(KT  + T ) = 


) - U(KI)  - u(K-l)l  , 


(126) 


This  expression  is  in  the  same  form  as  equation  (IO6)  discussed  earlier 


where : 


A = u(KT)  and, 

A°  = uttl 

1 1! 


and  where  u(t)  has  again  been  defined  as  in  equation  (104). 


(107) 


This  modified  scheme  of  control  is  mechanized  as  the  FOH  option  of 
the  program  in  Appendix  E and  is  siibstituted  for  the  ZOH  in  the  simu- 
lation by  specifying  FLAG  = 2 in  the  program. 


The  results  of  introducing  a First  Order  Hold  into  the  problem  are 
ejqjressed  by  two  typical  control  time  history  plots  idiich  follow.  Both 


plots  are  for  Q and  R weightings  of  unity.  The  first  plot.  Figure  33, 
is  a sample  rate  of  .01  seconds  while  Figure  34  is  for  the  sample  rate 
of  .02  seconds,  just  discussed  in  connection  with  the  ZOH  control  scheme. 

It  is  evident  from  Figure  33  that  the  control  has  been  smoothed 
considerably.  To  some  extent,  this  is  also  the  case  for  the  lower 
sample  rate  control  time  history  plot.  Figure  34.  In  this  figure, 
however,  the  "sawtooth"  character  of  the  FOH  can  be  seen  as  the  con- 
trol propagates  along  the  constant  slope  then  abruptly  jimips  when  the 
next  control  update  occirrs. 


Of  greater  significance,  however,  is  the  effect  of  the  FOH  on  the 
time  it  takes  the  system  to  reach  steady  state.  This  information  is 
presented  in  comparative  form  to  the  ZOH,  in  Table  VII. 


Table  VII 

ZOH  vs  FOH  Control  Scheme 

T = .02  R = Q = 1 

ZOH 

FOH 

lime  to  achieve  controlled 
steady  state  (sec) 

.556 

.630 

Maximum  positive 
deflection  (Rad) 

.005695 

.005054 

Maximum  negative 
deflection  (Rad) 

.01911 

.02718 

Final  control  value  at 
steady  state  (Rad) 

.001571 

.001571 
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It  is  apparent  from  Table  VII  that  although  both  hold  devices,  as 


mechanized  in  the  simulation,  achieve  the  same  steady  state  control 
value,  it  is  the  ZOH  which  achieves  steady  state  the  fastest.  Addition- 
ally, although  the  maximm  positive  deflections  of  the  control  surface 
are  comparable  (the  same  at  least  for  the  first  three  significant 
figiares),  the  maximiun  negative  excursion  of  the  FOH  exceeds  that  of 
the  ZOH  by  29.7  percent.  Though  the  same  results  are  achieved,  the 
FOH  is  not  as  energy  efficient  as  the  ZOH  in  light  of  the  greater  de- 
flection needed  to  control  the  aircraft  and  in  view  of  the  fact  that 
surface  deflections  produce  drag.  Longer  elapsed  times  to  achieve 
steady  state  and  greater  negative  excvirsions  of  the  control  were  the 
case  whenever  the  First  Order  and  Zero  Order  Holds  were  compared.  In 
view  of  this  sittiation,  the  ZOH  was  retained  as  the  most  efficient  of 
the  two  hold  schemes. 


Simulation  Observations 

Concentrating  on  the  ZOH  then.  Table  VIII  presents  those  combinations 
of  sample  rates  and  control  rate  penalty  weightings  (R)  investigated  while 
the  trajectory  error  weighting  (Q)  was  held  at  unity.  These  particular 
combinations  were  chosen  in  consideration  of  the  time  2ind  resources 
available  for  this  investigation, weighed  against  a desire  to  get 
reasonably  good  coverage  of  the  available  regime. 

Table  IX  presents  the  optimal  feedfoirward  gains  foimd  for  each 
case,  irtiile  Tables  X,  XI,  and  XII  present,  respectively,  the  three  op- 
timal feedback  gains  , and  also  detennined  as  optimal 

for  each  case  investigated.  Selected  points  from  Table  IX  and  X have 
also  been  presented  in  graphical  form  in  Figures  35  and  36,  respectively. 
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Table  VIII 

Svmnnary  of  Cases  Investigated 


From  Table  IX  and  its  associated  Figure  35^  it  is  seen  that  an  in- 
crease in  the  control  penalty  weighting  (R)  while  the  sample  rate  (T)  is 
held  constant,  produces  an  increase  in  the  value  (less  negative  value) 
of  the  gain.  However,  if  the  value  of  the  control  penalty  is  held 
constant  while  the  sample  rate  is  decreased,  the  optimal  value  of  this 
gain  decreases  (more  negative  value).  Similar  observations  can  be 
taken  from  Table  I with  its  associated  Figure  36  and  Tables  XI  and  HI. 
An  increase  in  the  penalty  weighting,  vrtiile  the  sample  rate  is  held 
constant,  has  the  same  effect  as  decreasing  the  sample  rate  while 
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Table  H 
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FIG  36/sAMPLE  RATE  vs  N,,^  (q-i,  var,  r) 


SAMPLE  RATE 


Table  HI 


4573  - .4556  - .4486  -.4434  - ./4364 


Table  XIII 


holding  the  control  penalty  weighting  constant.  The  result  is  an  ef- 


fective decrease  in  the  value  of  gain.  This  is  true  for  both  the  Nj 
and  feedback  gains.  From  Table  HI,  it  is  seen  that  an  increase 

in  the  control  penalty  weighting  (R),  while  the  sample  rate  is  held 
constant,  has  the  same  effect  as  decreasing  the  sample  rate,  while 
holding  the  control  penalty  weighting  constant.  The  effect  is  to  in- 
crease the  gain  . Finally,  Table  XIII  presents  the  time  required, 

6ii 

as  a function  of  sample  rate  and  control  penalty,  for  the  simulation  to 
achieve  steady  state  lock-on  to  the  commanded  response  with  zero  error. 
It  is  interesting  to  note  that  for  higher  sample  rates  (i.e.,  sample 
rates  greater  than  T = l/50)  there  is  little  change  in  this  amoxmt  of 
elapsed  time  required. 

These  observations  are  sxmmarized  in  Table  XIV.  Also  included  are 

•X 

the  observed  effects  of  these  operations  on  the  time  required  for  the  C 
system  output  to  match  C commanded,  with  zero  steady  state  error,  from 
Table  HII. 


Table  XIV 

A=  increase  Nummary  of  Gain  Effects  decrease 

Gain 

X.  Values 

Dperation  N. 

^d  \ 

Gadn 

Values  * 

Time 

RMuired 
for  C^  act= 
— C com 

Increase  R 
(T  Constant) 

A T T A 

▲ 

Decrease  T 
(R  Constaint) 

f T T A 

▲ 

* From  Table  XIII 
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Increasing  the  control  penalty  weighting  R then,  has  the  effect 
of  decreasing  the  size  of  the  controls.  This  has  the  same  effect  on  the 
system  as  increasing  the  sample  rate.  These  are  reasonable  effects  in 
that  increasing  the  penalty  on  using  excessive  control  rate  results  in 
smaller  controls,  while  decreasing  the  sample  rate,  or  the  frequency 
with  which  new  controls  are  introduced,  forces  the  system  to  adjust  by 
producing  larger  controls  in  order  to  keep  the  system  \mder  control. 

The  sample  rate  can  effectively  be  traded  off  against  the  control  pen- 
alty weighting  (R)  to  control  the  response.  The  time  required  to  reach 
steady  state  with  zero  error  increases  in  one  case  due  to  the  use  of 
smaller  controls.  In  the  other  instance,  with  decreasing  sample  rates, 
the  increase  in  time  is  due  to  the  associated  greater  drop  in  the  feed- 
back gains  Nj  and  Nj  verses  the  rise  in  gain  N,  and  the  use  of  larger 
values  of  ly. 


Saturation  Effect 

Observations  on  the  nature  of  the  controls,  especially  the  initial 
control,  indicate  that  a possible  problem  may  exist  due  to  the  large 
initial  jmp  that  occurs  in  the  control.  This  jump  in  the  control 
might  very  often  be  too  fast  and  lead  to  the  rate  saturation  of  a 
"real-world"  servo  actuator.  The  effect  of  the  20H,  as  can  be  seen 
from  Figure  31,  is  to  introduce  a series  of  steps  in  the  control. 
Consequently,  a series  of  impulses  would  be  characteristic  of  the 
control  rate  (u)  which  might  lead  to  rate  sat\iration  of  the  control 
actixator. 

The  maximum  movement  rate  of  the  horizontal  stabilizer  for  the 


TF-16  is  currently  60  degrees/seconds.  In  radian  measure,  this  is 
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1,047  radians/second.  Realizing,  of  course,  that  the  natxare  of  the 
simulation  is  a constraint  here  in  that  ,002  second  was  decided  upon 
for  its  iteration  rate,  it  is  nevertheless  valuable  to  characterize 
the  nature  of  this  possible  flaw  using  the  simulation. 

In  1/500^^  second,  a displacement  of  the  surface  of  ,00209 
radians  or  ,12  degrees  becomes  the  limiting  rate  of  movement  for  this 
control  surface.  This  limiting  rate  becomes  significant.  Consider 
the  case  in  vdiich  the  sample  rate  is  slowed  down  considerably.  During 
the  time  inteinral  that  the  aircraft  model  is  allowed  to  change  states 
between  updates,  the  plant  may  diverge  considerably  from  that  required 
to  null  the  error  between  C commanded  and  C actual,  Vftien  a new  con- 
trol is  calculated,  based  on  existing  states,  the  result  may  be  a con- 
trol correction  which  exceeds  the  capability  of  the  servo,  A.s  an 
example,  consider  the  case  in  which  the  present  control  is  -,0085 
radians.  The  control  update  rate  is  assumed  to  be  1/t  seconds.  After 
1/t  seconds  have  passed,  the  new  updated  control  is  calculated  to  be 
+,0065  radians.  The  absolute  range  of  control  can  be  seen  to  be  ,015 
radians.  The  servo  of  the  aircraft  is  in  effect  being  told  to  supply 
a ,015  radian  correction  in  l/500^  second.  Since  this  exceeds  the 
servo  limit,  this  control  cannot  be  supplied  indicating  that  the 
sample  rate  of  I/t  seconds  is  infeasible.  The  possibility  of  such 
an  effect  must  be  considered. 

For  this  simulation,  it  became  apparent  that  a limit  would  be  ap- 
proached in  terms  of  a minimum  sample  rate  ^diere  tiie  response  just 
falls  within  the  most  restrictive  of  the  C response  criteria  envelopes 
using  controls  vdiich  are  just  below  the  saturation  limit  of  the  servo. 
Such  a point  was  reached,  the  results  of  which  are  shown  in  Figure  37 > 
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38,  and  39.  The  minimum  sample  rate  was  found  to  be  T = 1/30  seconds 
with  a trajectory  error  penalty  weighting  of  unity  and  a control  penalty 
weighting  of  200.  The  system  C output  achieved  a steady  state  match-up 
with  the  commanded  response  in  1.716  seconds.  Greater  values  of  control 
penalty  or  slower  sample  rates  caused  the  response  to  be  stretched  out 
STifficiently  to  fall  outside  the  bounds  of  the  envelope  at  some  time  in 
its  time  history. 

Root  Migrations 

The  effects  of  variation  in  sample  rate  for  a constant  control 
penalty  and  the  effects  of  variation  in  the  control  penalty  for  con- 
stant sample  rates  are  presented  in  Fig\ires  40  and  41,  respectively. 

Both  figxu’es  portray  limited  data,  for  the  purpose  of  clarity,  in  the 
upper  right  (positive)  quadrant  of  the  Z-plane  unit  circle.  A con- 
jijgate  coimterpart  root  is  assimied.  Figure  40  depicts  the  concentric 
natiire  of  the  cxirves  of  constant  control  penalty  weighting;  the  inner- 
most curve  having  the  greater  weighting.  It  is  apparent  that  for  one 
particular  sample  instance,  such  as  T = l/lOO,  the  root  locations 
vary  considerably  depending  on  the  control  peneilty  weighting  selected. 

It  appears,  due  to  the  consistency  in  the  nat\are  of  the  overshoot  of 
the  output,  that  the  amount  of  overshoot  (i.e.,  f((i))  is  controlled 
by  the  trajectory  error  penalty.  With  Q fixed  at  imity,  the  system 
is  operating  in  a less  than  critically  damped  region,  in  the  neighbor- 
hood of  a .69  damping  ratio.  The  effect  of  a variation  in  the  control 
penalty  is  to  alter  the  particvilar  damped  natural  frequency  UJ^  . The 
direction  of  this  variation  in  natiiral  frequency  is  inversely  propor- 
tional to  the  direction  of  vauriation  in  the  control  penalty  weighting; 
that  is,  a decrease  in  03  results  from  an  increase  in  R.  This  appears 
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reasonable.  It  was  noted  earlier  that  increasing  the  control  penalty 
weighting  "stretched-out"  the  C*  response  of  the  system.  Such  an  ef- 
fect would  also  increase  the  time  required  for  the  response  to  reach 
its  peak  (tp).  Peak  time,  directly  related  in  this  way  to  the  control 
penalty,  is  known  to  be  inversely  proportional  to  the  natiiral  frequency 
thtis  completing  the  relationship. 

In  Figure  41,  target  roots,  based  upon  a i = .707  and  = .5.5 
as  presented  in  the  previous  chapter  have  been  calculated  and  super- 
imposed for  various  sample  rates.  The  locus  of  roots  for  two  parti- 
CTilar  sample  rates  are  shown.  Due  to  the  direct  mapping  nature  of 
Figure  27  to  Figure  28,  the  roots,  being  above  the  target  in  each 
case,  confirm  that  a damping  factor  less  than  .707  is  in  operation. 

C 
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VIII.  Conclusions  and  Recommendations 


Conclusions 

This  investigation  demonstrates  that  a discrete  optimal  controller 
for  the  longitudinal  pitch  axis  of  the  YF-16  Lightweight  Fighter  Proto- 
type  aircraft,  based  on  the  proposed  C transient  response  handling 
qualities  performance  criteria,  is  possible  at  .8  Mach.  A reduced 
state  model,  based  on  a short  period  approximation,  was  successfully 
developed  using  stability  derivatives  calculated  from  available  wind 
tunnel  information.  For  the  case  of  subsonic  flight  at  Mach  .8  at 
sea  level,  a positive  Cjjj^stability  derivative  was  shovm  to  exist. 

Such  a situation  is  indicative  of  a statically  imstable  aircraft. 

In  this  respect,  the  YF-16  is  unique.  Once  disturbed  from  equilibrium 
by  a pitch  up  command,  for  instance,  the  aircraft  will  continue  to 
pitch  up  (Fig.  12a). 

The  design  is  based  on  the  requirement  that  the  aircraft  success- 
fully -track,  with  zero  steady  state  error,  a 1-G  climb  step  input 
command.  Due  to  its  unstable  character,  results  confirm  that  the 
aircraft  attempts  to  continue  its  pitch  up  necessitating  the  use  of 
nose  dovm  elevator  control,  in  the  steady  state,  to  continue  the 
climbing  maneuver. 

In  answer  to  the  series  of  questions  vdiich  introduced  the  pre- 
vious chapter,  the  system  could  be  controlled  and  controlled  within 
the  innermost  C envelope  bounds.  The  everpresent  problem  of  rate 
saturation  was  shown  as  presenting  limitations  to  the  concept.  The 
tendency  for  initial,  large  jximps  in  control  position  was  shown  to 
lead  to  control  rates  which  could  exceed  the  movement  rate  of  the 
control  surface. 
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Through  the  use  of  a digital  simulation,  a ZOH  approach  was  found 
superior  to  a FOH  control  scheme  in  reducing  the  amount  of  elapsed 
time  for  the  system  to  achieve  a steady  controlled  state.  Investi- 
gation of  various  sample  rates  using  the  simulation  model  showed  the 
optima]  approach  to  a solution  only  slightly  effected,  in  terms  of 
elapsed  time  to  reach  steady  state  lock-on,  to  the  sample  rate  used. 
The  time  required  to  reach  steady  state  was  also  shown  to  be  appreci- 
ably imeffected  at  sample  rates  greater  than  T = 1/50  second.  More 
frequent  sampling,  however,  did  result  in  the  production  of  smaller 
magnitudes  of  control  deflections. 

Use  of  a trajectory  error  weighting  of  lanity  produced  an  over- 
shoot in  the  response  of  approximately  five  percent  associated  with 
a damping  ratio  of  .69. 

The  investigation  of  system  closed  loop  conjugate  root  migrations 
showed  the  trajectory  error  weighting  to  be  the  dominant  factor  in 
determining  the  amount  of  overshoot  and  the  effective  system  damping. 
Likewise,  the  control  penalty  weighting  was  seen  as  detennining  the 
system  natural  frequency.  This  relationship  held  that  the  natural 
frequency  was  inversely  proportional  to  the  control  penalty  weighting. 
Finally,  increasing  the  control  penalty  weighting,  while  the  sample 
rate  is  held  constant,  or  decreasing  the  sample  rate,  while  the  con- 
trol penalty  weighting  is  held  constant,  is  seen  as  causing  an  in- 
crease in  the  elapsed  time  to  steady  state  lock-on  of  system  response 
to  commanded  response. 

In  svnnmary,  digital  longitudinal  pitch  control  of  the  YF-16 
\ising  the  discrete  G*  approach  presented  here,  as  an  alternative  to 

the  present  analog  control  configuration  of  the  IF-16,  is  a viable 
possibility  warranting  further  investigation. 


115 


Recommendations 


r 


It  is  recommended  that  this  research  be  broadened  to  incorporate 
a larger  portion  of  the  YF-16  flight  envelope.  Stability  derivative 
information  for  three  additional  flight  conditions  is  presented  for 
futiire  research  purposes.  A full  scale  analysis  of  major  subsonic 
and  supersonic  operating  points  is  urged  to  confirm  the  applicability 
of  a discrete  C controller  and  the  nature  of  the  C response.  Quite 
possibly,  the  use  of  a system  of  modeling  stability  derivatives  based 
on  simplifying  assumptions  and  approximations  different  from  the 
Blakelock  approach  used  here,  could  be  useful  to  confirm  the  validity 
of  the  YF-lb  model  developed  in  this  investigation.  From  the  yiman 
variation  in  the  magnitude  of  the  controller  gains  with  sample  rate, 

I noted  in  this  study,  it  would  be  valuable  to  investigate  the  magnitude 

' of  the  gains  for  other  flight  conditions.  Relatively  close  gain 

, results  between  operating  regions  might  foster  the  identification  of 

1 

I one  fixed  set  of  gains  or  a limited  range  of  gain  values  which  con- 

i trol  the  system  over  a broad  range  of  operating  conditions.  Failure 

i 

to  identify  such  "all  purpose  gains"  might  necessitate  an  adaptive 
controlling  gain  scheme.  Additionally,  the  feasibility  of  substituting 
I the  optimal  gains  from  one  flight  condition  for  its  eqvdvalent  at  some 

1 other  operating  condition  should  be  investigated.  Such  a "mismatch 

i 

i analysis",  though  suboptimal  in  its  approach,  might  still  provide 

control  of  the  system  and  aid  in  an  all  purpose  gain  identification. 

This  study  concentrated  on  a step  input  pilot  command  to  the 
system.  Future  investigations  should  address  alternate  inputs  to 
the  system.  Additionally,  the  approach  taken  here  was  completely 
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deterministic.  Attention  should  be  given  to  the  discrete  stochastic 
aspects  of  a C**'  control  design  which  incorporates  random  distxirbances 
to  the  system. 


Finally,  a sensitivity  analysis  shovild  be  conducted  on  the 
effective  change  in  system  response  to  changes  in  the  values  of 
the  coefficients  , and  K"  in  the  original  C equation. 

Different  weights  to  these  terms,  effecting  the  blending  of  pitch 
rate  and  normal  acceleration,  might  yield  among  other  characteristics, 
a more  responsive  control  system. 
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Appendix  A 

Determination  of  Length  to  Tail  Jl^ 


This  appendix  first  defines  and  then  determines  the  value  of 
needed  in  the  evaluation  of  various  stability  derivatives  in  Chapter 
II.  Dimensional  values  are  taken  from  Reference  13. 

= distance  between  the  quarter  chord  point 
of  the  wing  mean  aerodynamic  chord  (MAC) 
and  the  quarter  chord  point  of  the  hori- 
zontal stabilizer  MAC. 

The  value  of  ji  was  determined  from  the  physical  dimensions  of  the 
aircraft. 


it 


jO 


Side  View 

IF-lb  Lightweight  Fighter  Prototype 


Appendix  B 

TRAMFUN  Pronrani  Inrut/Output 


The  determination  of  the  characterictic  equations,  characteristic 
roots  and  transfer  functions  from  the  three  longitudinal  equations  of 
motion  is  accomplished  using  the  TRANFUN  digital  computer  program 
(Ref.  16). 

For  Set  A,  the  equations  can  be  expressed  as: 


1.726s  + .0656  ’ .0005  .0621 

'1 
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0 

.1882  1.7325s  + 4.4942  -1.7073s 

-.4504 

0 .0093s  - .1788  .0135s^  + .0268s 

. 

M— 

-.6452 
L J 

This  information  is  input  into  TRAI-IFUN  in  the  following  manner: 
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The  following  results  are  forthcoming: 
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The  results  determine  the  system  transfer  fxmctions  which  can  be 
expressed  as: 


tC 


.00000304023^  .06971263s  + .18506976 

.040368983s^^  + .2137988s3  _ .3l093424s^ 


- .012018067s  - .0020896749 


-.25997  (s  + .0.1899276  + . 05968549 j)  (s  + 183.1493) 

(s  - 1.223932)  (s  + .02093674  ± .07804015j)  (s  + 6.478175) 


-47.6135  (s  + .03799498)  (s  + 2.676.139) 

(s  - 1.223932)  (s  + .02093674  ± .07804015j)  (s  + 6.473175) 


(131) 


(132) 


(133) 
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Phugoid  Roots  (Set  A.) 


The  phugoid  factor  of  the  characteristic  eqiiation  is: 

s + .02093674  ± .07804015J 

From  this  w < <■>  > 4 > and  T for  this  oscillatory  mode  can  be 
determined: 

% =V  «»/■  = *0808 

<«V  = .07804 

= .25912 

T:=  = 77.762  sec 

Short  Period  Root  (Set  A) 

The  short  period  factor  of  the  characteristic  equation  for  this 
unstable  situation  is: 

+ 5.254243s  - 7.928845684  = 0 
For  Set  B,  the  equations  can  likevd.se  be  expressed  as: 


’l.l506s  + .1033  .0813  .0276 

1 

0 

.1200  1.1465s  + 4.8132  -1.11862s 

<C 

T, 

- 

-.361 

• 

0 -.00164s  + .7392  .006s^  + .01855s_ 

-.5157 

This  information  is  input  into  TRANFUN  as: 


(134) 


(135) 


(136) 


(137) 
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The  follcv/ing  output  infoi'mation  results: 
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The  resulting  system  transfer  functions  are; 


(s  + .9612166)  (s  + 361.2223) 


.00791497743^  + .056298662s^  + 1.0590759s^  + .094453828s  + .0024482304 


•c 

X 


(138) 


-.3U87  (s  + .0/»48B49  ± .02300l64j)  (s  + 269.4224) 

(s  + .04474255  ± .01790868j)  (s  + 3.511721  + 10.99289j) 


(339) 
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(140) 


-86.03606  (s  + .08775803)  (s  + 3.745108) 

(s  + .0A474255  ± .01790868 j)  ( s + 3.511721  + 10. 99289 j) 


Phvigoid  Roots  (Set  B) 

The  phiigoid  factor  of  the  characteristic  eqiiation  is: 


s + .0A47255  ± .01790868j  = 0 


From  this  Ca  > Uy  , S , "17 

n • 

U)^  =V«-*+  ioj- 


= .0179 


Short  Period  Root  (Set  B) 

The  short  period  factor 


are  fomd  to  be: 

= .0482 

= .9279 
= 130.357 

the  characteristic  eqiiation  is: 


s + 3.511721  + 10. 99289 j = 0 


From  this  ^ > and  are  determined  to  be: 


= 11.540 


(141) 


(142) 


(143) 


10.99289 


Appendix  C 


Continuous  State  Variable  Equation 


This  appendix  details  the  development  of  a state  variable 
representation  of  the  longitudinal  equations  of  motion.  Such  a 
development  is  needed  in  order  to  model  the  TF-16  as  a linear 
system  of  the  form  x = A x + B u for  implementation  of  the 
controlled  system  model  represented  in  Figiire  19. 

Returning  to  the  original  equations  shown  previously  in 
Chapter  II: 


^ ^ ^ X ^ ^ 

-(L  i-a  u^2iUe._c_c  e.c^€>e=<l  r (2) 

^ Sg  c?2(  ^ ^ 5>.  «2>  > 

-c.,  i-lc  i.c.i  - 5 c © ^ (3) 

U Jj  .C-S  •>%  ✓ fc 


•*-  QtL  ^ ot 


jJt  / 


and  letting 


fe  A 

'W 

L A 


k.  = -2- 

^ Am 


^ C145) 


vrtiile  neglecting  the  following  terms  by  setting  them  equal  to  zero 
for  the  reasons  explained  in  Chapter  II: 


ft 


Cx.  = = 

Cyp  s <S>  r 0 

div 

(146) 

the  eqiiations  of  motion  can  now  be  written  as: 

. * / ✓ 
ka.c.  < 

1 VW 

(147) 

- - c,  i,k,o-k/i^e=  e> 

(148) 

w 

- ^3®"*  * - + fc*  * - ® ^ S, 

^Jv  = 

(149) 

Dropping  the  primes  and  defining  the  following  state  variables: 

— - 

Xt  = U 

X 

1 

= <X- 

X = 

= Q 

^ = q 

1 

(150) 

The  equations  become: 

k i - C.  X - X,  - 

. c y 

= 0 

(151) 

(152) 

j3 

1 - ^ 

{ 

(153) 

/» 

Dividing  equation  (151)  by  K,  yields: 

* £*»  X,  ♦ £*.<  x» 

- V 

•^i 

A =0 

(154) 
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1 


y _ y 4 y ^ v 

• I *^1 


(155) 


Similarly,  equation  (152)  can  be  rewritten  as: 


(156) 


or  after  expanding: 


- K - Kj  (157) 


Solving  for  X2 


k X 


(158) 


or  finally  as: 


N L J n 


Since  6 = q,  then: 


*3  = *4 


(I6l) 
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The  final  equation  results  by  dividing  equation  (153)  by 


K*.  '^4.  D Xjt 


Which,  when  solving  for  becomes: 


Substituting  in  the  previous  expression  for  ; 


(163) 


£!fi  y + .^bL  +/ iil^A 


^*^•1  y + — C^H,  Xy  — k fl 

A It  *"»  V ly  k 


(164) 


M\iltiplying  through  each  term: 


ejx.  +rAc?M(?Jx  +[A(ic+fcCi 


h.C^  c„1i  4 ^ + Jk.  f. 


0.65) 


or  finally: 


(166) 


1 


Svmmarizing,  the  foiir  state  variable  equations  corresponding  to  u, 

. • 

0(.  , Q,  and  0 are  respectively: 


X,  = 


y ^ ___2<  X-  -f-  


k ' 


— «■  y,  + 


I “•>  J S 


= 1.0  Xi. 


(155) 


(160) 


(161) 


x = djy,  +rj!ic  c.  + ^jy^  + 

r^/k  ♦ k Cj  \c^  + J^cLly,  (^h,  4 O £ 

r,  ^ is  introduced  as  the  new  state  x^,  since: 


(166) 


20 

S + 20 


(75  ) 


where  is  the  commanded  deflection  in  the  horizontal  stabilizer, 

then; 


Jo  * 20^1^ 


and  solving  for  ^ : 


(167) 


-20  + 20^^ 


(168) 


333 


or  finally  that: 


-joXj.  ♦ 


(169) 


77 


the  new  forcing  function  u(t)  become  «|,^(t). 

Now  the  system  modeling  equations  of  motion  in  the  form  of 
eqixation  (77  ) become: 


r 


'I 


It 


When  the  appropriate  terms  are  extracted  from  Table  III,  and 
substituted  into  this  latest  expression,  the  follovring  equation 
results: 


r 1 

* 

• • 

^1 

-.03800914  -.99928970 

-.03598122 

0 0 

*1 

0 

*2 

-.10863579  -2.5942135 

0 

.985425739  -.25998704 

0 

= 

0 0 

0 

1.0  0 

0 

.07495323  15.0520989 

0 

-2.6725778  -47.677365 

*4 

0 

*5 

0 0 

0 

0 -20.0 

*5 

20.0 

L.  J 

^ -J 

(171) 


which  is  positionally  equivalent  to  the  terns  of  x = A x + B u. 

Equation  ( 171)  can  be  simplified  in  terms  of  the  number  of  states 
by  using  a short  period  approximation.  Such  a simplification  is  justi- 
fied since  the  interest  here  is  in  the  transient  response  of  the  air- 
craft. Additionally,  the  coefficients  of  the  velocity  state  (u)  in  the 
A matrix  above,  have  only  a minor  contribution  on  the  short  period 
transient  response.  The  velocity  state  (u)  can  therefore  be  eliminated 
without  degrading  the  short  period  performance  of  the  aircraft. 

Again,  referring  to  the  short  period  approximation  equations 
of  motion  from  Chapter  III: 


Ms 

-2i2i  s - 

si* 

II 

(61) 

■ h 

J 
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[■  :k  ^ ‘]  ‘ ‘“’ 


and  eliminating  the  primes  and  s's,  eqixation  (6l)  becomes: 


^ i - Ci  t - 2^  e = 


solving  for  ^ 


4 /.o  S - £ 


or  eqiiivalently. 


(172) 


(173) 


i = —iL  ot  -f  i.o  e 4 £l^ 

k, 


(174). 


Likewise,  eqioation  (62)  becomes: 


- J_c.  i-c^  a - 

jr>t  * ^ 


Xc^  © 


(175) 


o(  + ‘^  © - k.  e 

-*  Jf  ^ **  •*  A 


(176) 
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and  solving  for  © becomes: 


• • 

e 


Ka.  ® Ka.  (177) 


Now,  if  the  expression  just  found  for  ^ is  siobstituted  into  this 
newest  expression,  6 becomes: 


The  remaining  state  equations  remain  the  same.  Specifically: 

© = q (179) 

and 

= -20.0  + 20.0  (168) 

In  state  variable  format,  equations  (l68),  (l74)>  (l78),  and  (179)  can 
be  expressed  as: 


(180) 


As  the  highlighted  area  of  this  matrix  equation  indicates,  a 
trivial  state  relationship  exists,  i.e,,  the  derivative  of  a state 
equals  its  derivative  and  this  state  makes  no  contribution  to  the 
remaining  state  expressions.  The  0 state  can  therefore  be  eliminated 
reducing  the  state  matrix  equation  even  further.  The  resulting  reduced 
short  period  approximation  state  variable  equation  vdiich  models  the 
YF-16  longitudinal  dynamics  is: 


0 


vrtiich,  when  the  appropriate  values  from  Table  III  for  M = ,8  at  sea 
level  are  substituted,  becomes: 


Ik 


Appendix  D 

Development  of  the  RecTxrsive  Opt j ma1  Control 

This  appendix  presents  the  development  of  the  expression  for 
the  recxirsive  discrete  optimal  control  u(KT)  used  to  control  the 
system.  The  development  is  a combination  of  the  approaches  pre- 
sented in  References  8,  12,  and  15. 

Development 

Given  the  discrete  system  of  equations: 


c 

3c(K  + 1)  T = A^i  X (KT)  + Bd  u(KT) 

(96) 

7 

y (KT)  = Cd  X (KT) 

(97) 

it  is  necessary  to  find  the  sequence  of  controls  that  minimize  the 
discrete  eqiiivalent  of  the  continuous  quadratic  functional  which 
follows; 

^ [^o  " u(t)^  R u(t)^  dt 

® (Ref.  12)  (183)  ■ 

where  is  an  arbitrary  step  input  (Z^  = C com).  This  continuous 
cost  functional  can  be  replaced  by  the  following  equivalent  discrete 
expression  where  and  AT  replace  ^ and  Jt  , respectively,  and 
u(t)  is  replaced  by  its  first  difference  equivalent.  Note  that  an 
infinite  time  cost  functionaO.  is  used  because  it  results  in  a constant 
gain  control  law.  The  discrete  cost  function  becomes: 
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^pK+l)T  - u(KT)j  %-AT^u(K+l)T  - u(KT)jjArj  ^ ^ 


or  equivalently  for  this  particular  problem: 


= E {(' 


*^com  ~ act; 


r ^ ( 


* r*  \ 

^ com  - act  ) + 


^u(K+l)T  - u(KT)j'^  Rjj^u(K+l)T  - u(KT)j|  (185) 

Let  Xq  be  the  solution  of  equation  (96),  that  is,  x(KT)  = x(K+1)T  in 
the  steady  state.  Then,  equation  (96)  becomes: 


K. 

x(K+l)T  = X KT  = 

Ad  x(KT)  + Bd  u(KT) 

(186) 

or. 

3E(KT)(r  - Ad) 

Bd  u(KT) 

(187) 

'9*7 

Now,  if  Xn  = lim  x(KT)  and 

k-*«* 

u = lim  u(KT),  equation 

k-*«» 

(187)  becomes: 

■ss 

*0  ^ ^ ~ ®d 

(188) 

or. 

^ - i]  ^ 

Bd  uo 

(189) 

I 


11 

i’T 

I-) 


r y 


w 


•I'T' 


>’6 


idiich  is  the  unique  solution  of: 


Xq  = Ad  Xq  + B^j  Uo 


This  implies  that: 


Xq  - Ad  Xq  = Bd  Uq 


(190) 


(191) 


Additionally,  using  equation  (97),  if  = y(KT),  then  by  equation 
(189): 


2o  = CdM  Ad 

= Cd  Xq 

Solving  for  u : 


[-[  ^ ]‘^  ®<i 


(192) 


“o 


(193) 


The  development  continues  by  defining  the  following  error  variables. 
Letting  the  transient  in  the  control  be  defined  as: 


v(KT)  = u(K+l)T  - u(KT) 


then  the  trauisient  is  u(KT)  is  defined  as: 


u(KT)  = u(KT)  - u„ 


which  Implies  that: 


(194) 


(195) 


Uq  = u(KT)  - u(KT) 


(196) 
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Finally,  the  transient  part  of  x(KT)  is  defined  as: 


x(KT)  = x(KT)  - 3^0 

This  allows  the  definition  of  *3^K+1)T  as  follows: 

x(K+1)t  = x(K+1)T  - 3co 

and,  u(K+1)T  as: 

u(K+l)T  = u(K+l)T  - Uq 
vdiich,  from  equation  (196),  becomes: 

'u(K+l)T  = u(K+l)T  -[u(KT)  - uJ^KT)] 

Regrouping  the  terms: 

^(K+l)!  = ITCKT)  + |^u(K+l)T  - u(KT^ 
and  using  eqioation  (194),  this  becomes: 

uiCK+Dl  = u(KT)  + ”v(KT) 

Finally,  let  v(K+1)t  be  defined  as: 

T(K+1)T  = u(K+l)T  - u(KT) 

Subtracting  3^  from  both  sides  of  the  identity  x(K+1)T  = 
x(K+1)T,  equation  (96)  can  be  written  as: 

5E(K+1)T  - 3!o  = x(K+l)T  - 

= A;i  x(KT)  + Id  u(KT)  - ^ 


(197) 


(198) 


(199) 


(200) 


(201) 


(202) 


(203) 


(204) 


Substituting  equation  (190)  for  the  on  the  right  side  of 

equation  (204): 

3^K+1)T  - Xq  = Ad  x(KT)  + Bd  u(KT)  - ^ Ad  + % u J 
= Ad  [x(KT)  - ^]  + [u(KT)  - uj 
but,  x(K+1)T  - 3^  = XCK+I)!  by  equation  (198),  or 


x(K+l)T 


[x(KT)  - Xo3  + Bd  [u(KT)  - uj 


(205) 


(206) 


Now  substituting  in  equations  (195)  and  (197)»  the  result  becomes; 


x(K+1)T 


Ad  x(KT)  + Bd  u(KT) 


(207) 


Equations  (202)  and  (207)  can  be  expressed  in  state  space  notation  as: 


x(K+1)T 


u(K+1)T 


Ad  Bd  x(KT) 


01  u(KT) 


Let  the  matrix 


^d  ®d  A.  I 

IX  _ _ = 

0 1 -L 


and  the  matrix 


v(KT) 


(208) 


then. 


x(K+1)T 

u(K+1)t 


^ ■u(KT) 


v(KT) 


Realizing  that  equation  (I84)  can  be  expressed  as: 


•^(u)^  = 


- y(KT)]  ^ Qd  [ Zo  - y(KT)]  + 
Ju(K+l)T  - u(KT)]  Rd  [ u(K+1)T  - u(KT)jj 


(209) 


(210) 


/ 

V*-. 


and  recalling  from  eqiiation  (97)  that: 


y(KT)  = x(KT) 


and  from  equation  (192)  that: 


^o  “ ^d  ^ 


then. 


( 97) 


(192) 


Zq  - y(KT)  = [ Xo  - x(KT)]  = C^j  [ -x(KT)]  (211 ) 

and  the  cost  functional,  equation  (184),  can  be  rewritten  as: 


•’U),  ' 


K = 0 


[u(K+l)T  - u(KT)J  ^ % [u(K+l)T  - u(KT^]  (212) 

Equation  (194)  can  then  be  substituted  into  equation  (212 ) with  the 


result: 


v(Kl)  ^ % aKT)j 

= (KT)  ^ Qjj  Cd  x(KT)  + v^(KT)  Rd  ^KT) 


K = 0 


(213) 


(214) 


^ T Cd*^  Qd  i^dk  x(KT) 

1 1?r(KT)  '5T(kt)]  T-_- 

^ -*  0 ,0  u(KT) 

= 0 I Jl.  . 


^(KT)  Rd  vtKT) 


C215) 
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F 


where , 


“d'''  Hi  “d  i ° 


I 


1 °J 


*^Ricatti 


(216) 


The  control  TCKT)  which  minimizes  this  cost  functional  is  of  the 
form: 


jn 


J>s 


5'^ 


t‘0 


\ • x(KT)  + K2  . u(KT) 

d d 


\h  h 1 • 

x(KT) 

L d -"dj 

u(KT)_ 

(217) 


where  and  K2  are  the  gains  obtained  using  the  positive  definite 
d d 

steady  state  solution  (P)  of  the  discrete  Ricatti  equation  in  the 
expression: 


Ki  1 K2 

d I d 


- -[r¥r  r'p# 


(218) 


The  discrete  Ricatti  equation  used  is  of  the  form: 

' Siioatti  + p $ - <§’■  pr[f^  pr+  Kj]  pTjfi 


(219) 


The  final  expression  for  the  control  law  is  obtained  by  manipulating 
equation  (217).  Assume  that  there  exists  a matrix  W such  that: 


>®d  = I 


(220) 


Using  equation  (207),  repeated  here: 
IdCK+DT  = Adx(KT)  + 8^11(0) 


(207) 
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it  follows  that: 


“ITCKT)  = x(K+l)T  - x(KT) 


(221) 


multiplying  through  by  VT  and  thereby  employing  equation  (220): 


VB^  'u(KT)  = W x(K+1)T  -V  x(KT) 


(222) 


"u(KT)  = W ^K+1)T  - W A , x(KT) 


(223) 


When  equation  (223)  is  substituted  into  eqioation  (217),  the  result  is: 


'v(KT)  = Ki^x(KT)+K2  W x(K+l)T  - M Ajj'x(KT)j 
= Fki  - K2^  W A^]  x(KT)  + ^2^  W x(K+1)T 


(224) 

(225) 


Adding  and  subtracting  K2  W x(KT)  to  the  right  side  of  eqxxation  (225), 

d 

results  in: 


T(kt)  = 1%  - K2  w A^ 

L d d J 


d d 


x(KT)  + K5  W x(KT)  + 
'^d 


K2  W x(K+l)T  - Kp  ¥ x(KT) 
d d 


(226) 


T(KT)  = [ - ^2  v(k^  - '^‘(KT)  + 

K2  ¥ [x(K+l)T  ~ x(KT)] 


(227) 


Equation  (227)  can  be  simplified  by  defining  such  that: 


-1^  Cd  = Ki  - K2  W (A^  - I) 

d d 


(228) 


14s 


I 

I 


ja'' 


1® 


ii 


J33 


Now  eqiiation  (227)  becomes: 

T(KT)  = -L^  x(KT)  + W x(K+l)T  - x(KT)] 

= -L^^C^[x(KT)  - ^jj+  W [?(K+1)T-  x(KT)] 


but  from  equations  (97)  and  (192): 


Cd  x(KT)  - Xq  = y(KT)  - Z 


then. 


v(KT)  = -L^  [ y(KT)  - Z^j  + W [ x(K+l)T-  x(KT)] 

Z^  - y(KT)  j + K2  W [ x(K+l)T-  x(KT)] 


= 


but. 


"3^(K+1)T-'?(KT)  = [x(K+l)T-  Xoj  -Jx(KT)-XQj 


= x(K+l)T  - x(KT) 


therefore. 


T(KT)  = Itl  I - y(KT)  j + K2^  W ^x(K+1)T  - x(KT)J 


Using  equation  (194),  this  can  be  rewritten  as: 


(229) 


(230) 


(231) 


(232) 


(233) 


u(K+l)T  - u(KT)  = ly  f*  Z - y(KT)']  + K2  W [x(K+l)T  - x(KtJ 

o U ° J d **  /, 


(234) 


149 


Now  letting  Ko  • ¥ = N , , this  becomes: 


u(K+l)T  - u(KT)  = [ Zq  - + % [ ^(K+1)T  - x(KT^ 

where  Zq  is  the  system  input  , and  y(KT)  by  equation  (97)  is 

C x(KT).  Since  from  equation  (228): 


(236) 


equation  (220 ) can  be  expressed  as: 


-1,  . 


Using  this  information,  can  be  expressed  as: 

% = \\  - \ (*d  - i)  5d]  [^d  ®d] 


(237) 


(238) 


while , 


Nd  = K2 


_ V'  ( \ " “-o  ^ '] 


(239) 


"d  - {\  ^ ^d  ^d)  (*d  - -) 


(240) 


Eq\xation  (235)  is  the  resulting  recursive  control  expression  which,  when 
incorporatea  with  equations  (238)  and  (24o)>  produces  the  sequence  of 
optimal  controls  for  the  system  of  equations: 

x(K+l7T=  Ad  x(KT)  + B^j  u(KT) 
y(KT)  = IJd  x(KT) 
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(96) 

(97) 


Appendix  E 

Simvilation  Computer  Program 


The  foUovri-ng  program  (SIK3)  was  developed  and  used  for  this  in- 
vestigation  to  solve  the  C discrete,  optimal  control  problem  discussed 
in  the  text.  The  reader  is  referred  to  Figures  22,  24,  and  25  of  the 
text  for  a detailed  representation  of  the  flow  of  logic  in  this  algo- 

) 

rithm.  Explanatory  comments  are  included  in  the  program  to  clarify  key  1 

areas/options.  The  program  is  easily  adaptable  to  other  aircraft  or 

other  flight  conditions.  As  presented,  the  continuous  aircraft  plant 

(3  X 3)  dynamics  matrix  (A)  and  (3  x l)  control  matrix  (B),  as  well  as  I 

the  (l  X 3)  observer  matrix  (CD)  are  for  a YF-lb  flying  at  Mach  .8  at  ! 

sea  level.  Additionally,  the  Q and  R weightings  are  set  to  unity. 

The  program  was  run  on  the  Intercom-Scope  time  sharing  access  i 

system  at  AFIT.  For  execution,  the  program  requires  library  access  j 

I 

to  the  "AFITSUBROUTINES"  package  of  computer  programs  and  the  Aero-  1 

medical  Research  Laboratory  "CONTROL,  CY=4."  computer  subroutine  ■ 

package  both  of  which  are  available  at  AFIT.  j 
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Appendix  F 

Zero  Order  Hold  Simulation  Output 

The  abbreviated  computer  product  which  follows  is  the  output  of 
the  program  found  in  Appendix  E with  the  ZOH  option  selected  (FLAG  = 1). 
This  particular  run  is  for  a specified  sample  rate  of  T = l/AOth  second 
with  R = 150.0  and  Q = 1.0. 

The  first  page  of  the  output  indicates  that  the  sample  rate  has 
been  adjusted,  as  per  Table  VI,  from  the  specified  value  (.025)  to  the 
adjusted  value  (.026)  used  in  the  actual  simulation.  This  is  evident 
on  page  170  where  it  is  apparent  that  the  first  sample  occurs  after 
.026  seconds  have  elapsed.  The  control  applied  (u  = -.1819E-02)  is 
retained  until  the  next  sample  occurs  at  .052  seconds  at  which  time 
a new  control  is  calculated  and  applied  on  the  next  iteration. 

The  three  system  states  a , 0 , and  6^  are  listed  for  each 
time  increment.  Additionally,  the  system  output  (C*)  is  listed  along 
with  the  pilot  commanded  C response.  By  comparing  the  actual  C re- 
sponse to  the  C*  commanded  value,  the  tracking  capability  of  the 
particular  system  configuration  can  be  determined.  Successful  system 
"lock-on"  to  the  commanded  1-G  climb  is  evident  on  page  17I.  At  an 
elapsed  time  of  I.6O6  seconds  from  the  time  the  climb  command  (C* 
command)  is  initiated,  zero  error  is  achieved  and  maintained  by  the 
system. 

The  abbreviated  example  of  the  simulation  output  is  concluded  on 
page  172  vdiere  the  Z-plane  system  roots  are  listed. 
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THE  ST^IHTrOH  •RESULTS  FOLLOW 
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I Appendix  G | 

I Plotting  Algorithm  | 

! 

» J 

j The  following  program  (GRAPH)  was  developed  and  used  to  produce  I 

i the  plotted  information  fovmd  in  this  investigation.  All  plotted  1 

i 

information  was  produced  on  a ZETA  Research  Inc.  230  series  Zeta  | 

plotter  using  ZETA  PLOT  subroutines  (Ref.  19).  The  information  is  | 

supplied  by  the  main  program  found  in  Appendix  E.  The  main  program,  i 

in  addition  to  listing  the  information,  places  all  the  generated  data  j 

j 

onto  magnetic  tapes.  TapeS  is  used  to  store  the  values  of  the  three  I 

states,  C i ^ Com  > particular  control,  for  each  l/500th  ■ 

second  time  increment  of  the  simulation.  Tape9  is  used  to  store  an 
integer  which  equals  the  total,  nianber  of  lines  of  information  to  be 
found  on  TapeS.  For  example,  for  a one  second  simulation  run,  TapeS 
would  contain  501  lines  of  output  information  (includes  zero-time 
initial  condition  values)  and  Tape9  would  store  the  number  501.  As 
presented,  the  routine  is  limited  to  a meorlmum  of  1,001  lines  of  data 
corresponding  to  a maximum  simulation  run  time  of  2.0  seconds.  This 
limit  is  a result  of  array  and  field  length  restrictions  on  the  Inter- 
com-Scope time  sharing  access  system  on  which  the  program  was  run. 

; The  program  uses  the  subroutines  ZTPLOT  and  CCAU]^  both  available 

: at  AFIT,  which  must  be  attached  before  the  program  can  be  executed. 

i 

I On  the  intercom  system,  these  sxjbroutines  are  attached  as  follows: 

k 

i ATTACH,  CCAUI,  ID  = X65A321,  SN  = AFIT 

ATTACH,  ZTPLOT,  ID  = AFIT. 

LIBRARI,  ZTPLOT,  CCAUI 


j, 


w 


Additionally,  the  program  includes  a flag  which  detennines  which 
plots  are  produced.  This  is  accomplished  by  setting  the  IFIAG  variable 
in  the  program  to  1,  2,  or  3.  These  options  have  the  following  effects: 


Option 

Plots  Produced 

IFLAG  = 3 

States  vs  Time 
Output  vs  Time 
Control  vs  Time 

IFLAG  = 2 

Output  vs  Time 
Control  vs  Time 

IFLAG  = 1 

Control  vs  Time 

Nianerous  explanatory  comments  are  included  throughout  the  complete 
listing  of  the  program  which  follows: 
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